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Abstract 

In order to improve the understanding of particle vitiation effects in hypersonic propulsion test 
facilities, a quasi-one dimensional numerical tool was developed to efficiently model reacting particle-gas 
flows over a wide range of conditions. Features of this code include gas-phase finite-rate kinetics, a global 
porous-particle combustion model, mass, momentum and energy interactions between phases, and 
subsonic and supersonic particle drag and heat transfer models. The basic capabilities of this tool were 
validated against available data or other validated codes. 

To demonstrate the capabilities of the code, and to provide initial insight into the effects of various 
particle laden flows on ignition, a series of computations were performed for a model hypersonic 
propulsion test facility and scramjet. Parameters studied were simulated flight Mach number (Mach 5, 6, 
and 7), particle size (10, 100, and 1000 pm diameters), particle mass fraction (single particle and 1 
percent) and particle material (alumina and graphite). For the alumina particles, it was found that the 
presence of particles up to 1 percent mass fraction had very little effect on the gas phase, even though 
only the 10 pm particles closely followed the gas flow velocity and temperature. With the graphite 
particles, the 1 0 pm particles were either quickly quenched, or were quickly consumed, depending on the 
gas temperature. As the particle size was increased to 100 pm, the particles did not quench, but were still 
typically consumed within the model test facility. For the 1000 pm particles, combustion was diffusion 
limited, so particle and gas temperature had little effect on the combustion rate. When the particle mass 
fraction was increased to 1 percent, the main change was the addition of significant heat release. In those 
cases where low graphite reaction rates were observed for single particles, the increase to 1 percent mass 
fraction had very little impact. 

Hydrogen/vitiated air ignition delay calculations for the 1 percent mass fraction of graphite particles 
cases showed significant decreases in ignition delay in cases where higher graphite combustion rates were 
observed. Further calculations showed that this was due primarily to increased combustor inlet 
temperature, not the gaseous or solid vitiate species present in the flow. 

1.0 Introduction 

The principle challenge associated with ground testing of hypersonic airbreathing propulsion systems 
is the requirement to heat the air supply to temperatures in excess of 1300 K prior to expansion through 
the facility nozzle. Over the last 40 years, several different methods have been employed to accomplish 
this heat addition. The most commonly used method is combustion heating, wherein fuel is burned in the 
air stream to bring the flow temperature up to the required stagnation temperature, and then makeup 
oxygen is added to bring the flow composition back to the normal 2 1 percent oxygen mole fraction. 
Several facilities have been built which utilize an electric arc to heat the air in a stilling chamber to very 
high temperatures. This flow is then mixed with additional high pressure air to create an air flow at the 
desired stagnation temperature which is then supplied to the facility nozzle. Lastly, heat sink storage 
heaters have also been built for several additional facilities. In these facilities, a large mass of a refractory 
material in the form of pebbles or drilled blocks is heated to a temperature significantly above the desired 
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air stagnation temperature prior to the propulsion system test. The test air then flows through the 
refractory material which heats the air to a temperature slightly above the required stagnation 
temperature. Additional high pressure air is then added to the heated air to achieve the test stagnation 
temperature set point. It should be noted that heat exchangers and electrical resistance heaters that are 
typically used for lower Mach number propulsion test facilities cannot be used at the significantly higher 
temperatures required for hypersonic propulsion testing due to structural and heat flux limitations. 
Detailed descriptions of the operational procedures and limitations of each type of hypersonic propulsion 
test facility can be found in the JANNAF Scramjet Testing Recommended Practices and Guidelines 
(Weber, 2002). 

Each of the heating methods described above introduces contamination into the facility air flow, 
typically referred to as vitiation. This flow vitiation results in some level of departure of the ground test 
results from that which would be seen in free flight. For the combustion heated facility, the products of 
combustion remain in the air flow. If hydrogen is used as the heater fuel, the product is primarily water. If 
a hydrocarbon fuel is used, carbon dioxide would also be present. If the fuel/air mixing scheme is not 
adequate or if the heater combustion length is insufficient for the given test condition, unbumed fuel may 
also be present in the flow. As the simulated Mach number is increased, the higher temperature of the 
mixture results in an increasing mole fraction of combustion radicals, such as H, O, and OH, as well as 
NO and N0 2 , exiting the heater and entering the nozzle, in addition to the increasing fraction of water 
vapor and carbon dioxide. If a hydrocarbon fuel is used, the fraction of carbon monoxide present in the air 
stream also increases with increasing simulated Mach number. The rapid expansion of the heater outflow 
through the facility nozzle freezes the composition at a point downstream of the nozzle throat, so that the 
mixture is not allowed to equilibrate. Under some circumstances, water droplets may form as the water 
vapor becomes supercooled due to the expansion to high altitude flight conditions. This condensation is 
rarely seen at the test point, in practice. 

As the combustion heated air enters the engine inlet, the changes in the ratio of specific heats, y, and 
the gas molecular weight alter the shock structure from what would be seen in clean air. A variety of other 
physical characteristics of the flow are also altered due to the flow vitiation, such as heat transfer rates, 
vibrational energy relaxation rates, turbulence, and inlet flow distortion, among others. As the flow enters 
the combustor and fuel is added, the combustion kinetics can be significantly altered, effecting both 
ignition and flame holding in ways that can be hard to predict, due to a variety of competing mechanisms. 
Lastly, nozzle expansion is also effected by the flow composition. As can be seen by the list above, 
combustion heater vitiation adds a large amount of uncertainty to the test results. Additional detail on 
these effects can be found in the summary papers of Powell and Stallings (1998), Pellett, et al. (2002) and 
Fry (2004). Additional difficulty is incurred with combustion heated facilities in terms of determining 
which primary variables should be matched during the testing. The changes in gas properties makes it 
impossible to match both total temperature and total enthalpy, both Mach number and velocity, or both 
total pressure and static pressure, for instance. The experimentalist must evaluate which properties are 
key(s) for each test, a subject of much debate in the hypersonic engine testing community through the 
years. 

If an arc heater is used, high levels of NOx are formed in the vicinity of the arc, which can persist in 
the mixed flow entering the nozzle at mole fractions above 3 percent. These compounds have little effect 
on the shock structure or other aerodynamic features, but can significantly catalyze the ignition process 
(Slack and Grillo, 1977). 

Lastly, a storage heater generally represents the lowest level of chemical vitiation. Like the arc heater, 
a storage heater can generate some NOx, but not to the same level. At most, an equilibrium concentration 
of NOx might be generated, depending on the residence time of the air mixture at high temperature. As 
the thermal generation of NOx is fairly slow, it is possible that NOx levels below equilibrium might be 
present in some facilities, depending on the specifics of the facility configuration. 

The above discussion gives a brief overview of the gas-phase vitiation issues for hypersonic 
airbreathing propulsion ground test facilities. However, each of these facility types can produce 
particulate vitiation as well. For the combustion heated facility, the particulate source is most likely soot 


N ASA/TM— 20 10-216765 


2 



formed during the combustion process due to the existence of fuel rich regions in the hydrocarbon fuel/air 
mixing zone. Additional particulates may be present due to piping corrosion caused by hot water vapor. In 
an arc heated facility, electrode erosion releases a small amount of metallic particulates into the air 
stream. In a storage heater facility, eroded particles from the refractory material, a ceramic or graphite, are 
found in the air stream. While the quantities and size distributions of these particulates have generally not 
been characterized, it can be assumed that the level of particulate vitiation is highest with storage heaters 
due to the extremely large surface area of non-surface abrasion resistant material exposed to the flow. 
While the most recent surveys of the state of the art regarding the effects of vitiation (Fry, 2004; Pellett, et 
al., 2002) acknowledge the importance of particulate vitiation, only a single study has ever been 
conducted of the potential effects of particulates (Mitani, 1995). In contrast, the surveys of Fry (2004) and 
of Pellett, et al. (2002), reference dozens of past studies addressing various aspects of gas-phase vitiation, 
reflected in the summary of effects described above. 

Mitani (1995) performed a series of theoretical analyses to determine the effect of low micron-sized 
inert particles, assumed to be in thermal equilibrium with the carrier gas, on scramjet combustor ignition 
and flameholding. Mitani concluded that a mass fraction above 0.1 percent can inhibit ignition through 
the mechanism of radical termination on the particle surface. This mechanism is only valid for particles 
smaller than 3 pm in diameter. For the particles to have a thermal heat sink effect, the study concluded 
that a mass fraction in excess of 1 0 percent would be required. 

In order to look at the effects of larger particles that are not in equilibrium with the carrier gas, as well 
as reacting particles, throughout the test facility flow path, a numerical analysis is necessary. Flowever, up 
to this time, no numerical modeling has ever been done to examine the potential effects of particle 
vitiation in hypersonic propulsion ground test facilities. Reasons for this include: 

1 . The effects of particles have been conveniently assumed to be of secondary importance compared to 
the effects of gas phase vitiation. 

2. The numerical analysis tools typically used for hypersonics flows do not contain particle flow models. 
The existing tools are multi-dimensional codes with finite-rate chemistry capability characterized by 
long run times, and so the addition of particle flow models is not attractive as it would tend to require 
even longer computational turn-around times. 

3. Particle-gas flow interaction models that cover the broad range of conditions experienced in a 
hypersonic propulsion test facility are not readily available. 

4. The mass loading and size distribution of particle vitiates has not been determined, so key inputs for 
model formulation are not available. 

While the first assumption, that particle vitiation effects are of secondary importance compared to 
gas-phase vitiation effects, may well be valid for combustion heated facilities, and perhaps arc -heated 
facilities as well, it is not a good assumption for storage heater-based facilities. In storage heater-based 
facilities, the amount of particle vitiation is likely to be substantially higher than in other types of 
facilities, and the amount of gas-phase vitiation is significantly lower. The final three reasons given above 
can be addressed, all or in part, through the development of appropriate particle-gas flow models and 
applying them to a quick turn-around numerical tool. By integrating the typically computationally 
demanding particle flow analysis into a more rapid analysis tool, parametric studies of particles of various 
sizes and other pertinent characteristics can be performed in a timely manner to understand under what 
conditions particle vitiation may be of concern. 

In this work, a quasi-one-dimensional (Q1D) numerical modeling tool and associated sub-models are 
developed that allow for a relatively rapid exploration of the effects of particle vitiation throughout the 
operating envelope of a storage-heated hypersonic propulsion test facility. A wide range of particle sizes 
are considered for both chemically inert particles (alumina) and reacting particles (graphite). In order to 
appropriately define the ranges of the parameters to be modeled, both a test facility and a scramjet engine 
model (up through the combustor) are laid out as part of this study. The test facility model is a 
compromise between the geometry and flow characteristics of the two large scale storage heater scramjet 
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test facilities found in the United States. This same model facility and model scramjet will then be used to 
demonstrate the use of the Q1D code for simulating hypersonic propulsion test flows with both inert and 
reacting particles. 

The first storage heater based hypersonic propulsion test facility is the GASL Leg IV facility. The 
facility consists of an alumina matrix storage heater followed by hydrogen fueled booster heater. The 
storage heater can heat dry air to a total temperature of 1390 K, simulating approximately Mach 5, with 
the booster heater then adding thermal energy up to a Mach 7 condition (Roffe, et al., 1997). The basic 
configuration is shown in Figure 1. 

The second configuration to be studied is the NASA Glenn Research Center Hypersonic Tunnel 
Facility (HTF) located at the NASA Plumbrook Station in Sandusky, Ohio. This facility consists of a 
drilled-core graphite block heater capable of heating gaseous nitrogen to 2500 K. The heated nitrogen is 
then mixed with ambient temperature gaseous oxygen to provide simulated air to the facility nozzle at a 
total temperature up to 2170 K (Perkins, et ah, 1998). This basic configuration is shown in Figure 2. 


Air 

Shutoff 

Valve' 


Air Supply 


Diluent 
Air Supply 


Booster 

Heater 


a 


/ 


Settling 

Chamber 


Pebble Bed 
Heater 


Facility 

T 


Figure 1. — GASL Leg IV diagram. 
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While the GASL Leg IV facility is vitiated above Mach 5, the model facility will assume non- vitiated 
flow for both the inert and the reacting particles at all test Mach numbers, as vitiated flow for the inert 
particle case would significantly affect the nozzle flow and inlet shock structure, making comparisons 
difficult between the different particle types. 

Neither GASL Leg IV nor HTF have ever been characterized as far as particle size and concentration 
as a function of test condition. Therefore, these particle properties will be dealt with parametrically when 
the code is exercised through a range of test conditions. There is no reason to assume that particle size is 
limited to a few microns diameter, as assumed by Mitani (1995). In the case of the HTF facility, particles 
on the order of 1 mm have been observed adhering to the heater graphite blocks during facility 
maintenance. Therefore, particle diameters in the range of 10 pm to 1 mm will be considered. Neither 
inert nor reacting particles with diameters less than 10 pm were used in this analysis as such small 
particles would be expected to track the gas flow very closely (Maxwell and Seasholtz, 1974) and, as will 
be shown in the following analysis, the 10 pm graphite particles almost immediately quench or are 
consumed, so a smaller particle would not be expected to show any different behavior. A particle mass 
fraction of 1 percent of the total test flow will be used as an upper bound, since a higher mass fraction 
would represent such a significant mass loss from the storage heater that heater refurbishment would 
likely be required after each test program. Such heater degradation has never been reported. While the 
results of this parametric study will not conclusively determine whether or not particle vitiation is of 
concern at any particular test point for a particular test facility, it will provide insight into what particle 
concentrations and sizes may be of concern for future research test programs. 
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In the following sections, the Q1D code and its sub-models will be described and the results of a 
series of validation cases based on the above range of conditions will be given. The code is then exercised 
through a parametric study of particle and gas flow combinations appropriate to the model test facility and 
model scramjet. Results are presented and conclusions are drawn as to the potential effects of particle 
vitiates. Finally, recommendations for further study are presented. 

2.0 Code Description 

2.1 General Approach 

The need to model a large number of facility and test hardware components over a wide range of 
parameter space, as well as the need to consider two phase flow with mass, momentum and energy 
interactions between the phases and finite rate combustion kinetics, limits the number of dimensions that 
can be practically used for the analysis. For this reason, a quasi-one-dimensional (Q1D) analysis approach 
was selected for the underlying gas-phase code. No existing code was available with the requisite 
characteristics for this study. A custom Eulerian (fixed grid), gas-phase, multispecies, time -accurate, Q1D 
computational fluid dynamics (CFD) code with finite rate kinetics was therefore written and validated. To 
this code was then added a Lagrangian particle tracking scheme with a global particle reaction model. 
Since larger particles are considered, an internal one-dimensional grid is generated for each particle being 
tracked which is solved at each time step for the particle internal temperature gradient and the heat 
transfer to the gas phase. If the particles are sufficiently small such that the Biot number is expected to 
remain less than 0. 1 (based on the particle radius as the characteristic dimension, (Incropera and Dewitt, 
2002) throughout the simulation, the particle is treated as isothermal. This is easily accomplished within 
the code by setting the number of internal particle grid points to two, as the number of shells within the 
particle is equal to one less than the number of internal grid points. Mass, momentum, energy and species 
source terms are included in the gas-phase governing equations to provide feedback from the particle 
drag, heat transfer, and reaction models. For inert particles, the particle reaction model is simply turned 
off. When the particle mass fraction is high enough that the number of particles being tracked would be 
too high for efficient and timely computation, each particle tracked is defined as being a representative 
particle for a larger number of particles at that same location. The number of particles represented by each 
tracked particle is provided as a user input. The mass, momentum and energy contributions from each 
tracked particle are then multiplied by the number of particles in each set to get the gas-phase source 
terms. 

This modeling approach allows for a large number of computations to be performed, but with certain 
drawbacks. Boundary layer related effects cannot be captured, nor can oblique shocks. Perhaps the 
biggest drawback for this problem is that the lateral diffusion of the particles when passing through a 
divergent section at high flow velocity is not captured. Previous studies have shown that when two-phase 
flow passes through a convergent-divergent nozzle, the particle stream will “separate” from the wall, 
creating a particle-free zone near the wall (Chang, 1980; Ishii, et al., 1987). The concentration of particles 
therefore increases in the rest of the flow. This effect is more pronounced as particle size is increased. 
However, since this study is parametric in nature, this effect can be accounted for when the results are 
applied to a specific configuration by correcting the particle concentration based on the nozzle flow 
conditions and the particle size under consideration. 

2.2 Gas Phase Governing Equations 

The gas phase governing equations were adapted from Kuo (1986), re-derived to include area 
variation. Viscous dissipation and diffusion terms determined to be unimportant for these simulations are 
omitted. The <|) term represents the local volume fraction of gas in the gas-solid mixture. The source 
terms, S m , S p , S a come from the solid phase relations, while the species source terms, S mJ , are a 
summation of both the gas phase and solid phase reactions. 
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Continuity: 
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Species Continuity: 
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where 
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Vi = 


A ^ 

Yi dx 


(7) 


The q w wall heat transfer term is currently a user input that can be used to “tune” the computational 
results to match available experimental data. If data are unavailable, an approach based on the Reynolds- 
Colbum analogy relating wall friction to heat transfer has been shown to be effective for Q1D 
calculations (Paxson, 1993) and could be readily incoiporated. 


2.3 Thermodynamic Properties 

The basic thermodynamic properties for each species, enthalpy (h g j), enthalpy of formation (h ro j), 
specific heat (c p J), Gibbs energy of formation (g/f) and molecular weight (Mi), were all taken from the 
JANNAF Thermochemical Tables (Chase, et al., 1986). The temperature dependent properties were 
piece-wise curve fit over four temperature ranges using third-order polynomials. The four temperature 
ranges were 0 to 300 K, 300 to 600 K, 600 to 1100 K, and 1 100 to 3000 K. The use of four temperature 
ranges provided high accuracy over the entire temperature range, typically within 0.1 percent of the 
original JANNAF data. From this data, along with the mass fractions of each species, the mixture 
properties can be defined, as follows. 
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For some of the calculations, the mole fractions of the chemical species are needed instead of mass 
fractions, which can be easily calculated as follows. 



Once the mole fractions are calculated, the mixture molecular weight can be obtained. 


(11) 


ns 

M = Y, X i M i 

i = 1 


This leads to the simple relation for the mixture gas constant, given the universal gas constant, 


R = ^-, 

M 


and the ratio of specific heats can then defined as 


Y = 


C p -R ' 


( 12 ) 


(13) 


(14) 


2.4 Transport Properties 

Species effective collision diameter (a,) and effective temperature ( T e i ) data, taken from Svehla 
(1962) and Hirschfelder, Curtiss, and Bird (1964), were combined with the Prandtl number (0.706) to 
generate the transport properties used in the calculations. Following the methodology found in White 
(1991), the binary diffusion coefficients are first calculated. 
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Where the collision integral is defined as 
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Once the binary diffusion coefficients have been calculated, an “average” diffusion coefficient for 
each species diffusing into the rest of the mixture can be calculated, using the following relation from 
Wamatz, et al. (1999). 


A = 


1 -Y, 


X: 


I Z). 

j=hj*i l ’J 


( 20 ) 


While the use of an average diffusion coefficient is somewhat less accurate than using a full multispecies 
diffusion model, it is much less computationally demanding and was deemed sufficiently accurate (~10 
percent) for the purposes of this study. 

To calculate the individual species viscosity, the method described in Hirschfelder, Curtiss, and Bird 
(1964) is adopted. 


Hi = 


26.993 (MiTgfi 


^ 2 , 2 &i 

where the collision integral, Q 2 , 2 > is represented by the curve fit 


( 21 ) 


Q 2 ,2 = 0.5015 + 2.5715(e“ zl448r * 2 .2 )+ 0.639 l(e“ a3856r * 2 .2 )+ 0.3546(e“ 0 016ir * 2 .2 ) (22) 


and 


T* 2,2 =~ 


The mixture viscosity is then calculated using Wilke’s law, as given by Diummond (1990). 
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Lastly, the mixture thermal conductivity is simply calculated by the standard expression 

k _ C p \i 
g Pr 


All of the above expressions assume laminar flow. 


(25) 


(26) 


2.5 Finite-Rate Kinetics 

Each of the N chemical reactions used in the gas-phase reaction mechanism, shown symbolically in 
Equation (27), is modeled using the standard Arrhenius expression for the forward reaction rate, 
Equation (28). 


N , N „ 

X v ' S ‘ ( 2 y ) 

/= 1 i = 1 


k f = AfT g n e R “ T g (28) 

The reverse reaction rate can be found from the relation 

kf 

h=~^~ (29) 

eq 

where the equilibrium constant, K eq , is calculated from the Gibbs energy of formation for each species and 
the reaction coefficients from Equation (27), following the procedure found in Kuo (1986). 

For this study, it was necessary to assemble a chemical kinetic mechanism that included H, O, N, and 
C compounds in order to capture the basic hydrogen combustion process for a scramjet combustor as well 
as the effects of the flow vitiates of interest. Based on a review of the available literature on gas-phase 
vitiation, as outlined in the introduction, the H-0 species H 2 , 0 2 , H 2 0, OH, H, O, H0 2 , and H 2 0 2 were 
selected. The nitrogen species N, HNO, NO and N0 2 were added to the basic diluent N 2 to capture 
dissociation and vitiation effects. Lastly, CO, C0 2 , and HCO were included as gas-phase products from 
the combustion of the carbon particles under study. The mechanism developed by Davis, et al. (2005), for 
H 2 -CO combustion was adopted for the bulk of the mechanism used for this study. The nitrogen 
chemistry was added from the GR1 3.0 mechanism (Smith, et al. 2005). The details of the combined 
mechanism are given in Appendix B. 

Once obtained, the forward and reverse reaction rate constants are used in the standard law of mass 
action expression to obtain the change in moles of each species for each reaction. These molar changes 
are then summed across all of the reactions, and then utilized as source terms within the governing 
equations. 
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2.6 Solid Phase Governing Equations and Numerical Solution Method 

There are two fundamental approaches that can be taken when formulating the governing equations 
for the solid phase in a two-phase computation. The first is to assume that the particles form a continuum 
that can be treated in the same manner as the gas phase. A parallel set of conservation equations to the gas 
phase equations can then be set up for the solid phase, simply exchanging (l-i)>) for <j) as the volume 
fraction occupied by the flow, and making the source terms of opposite sign in Equations (1) through (3). 
This is the approach, referred to as Eulerian, taken by Kuo (1986), Ludwig and Roth (1997), and Di 
Giacinto et al. (1982), among many others, and is typical of computations of fluidized beds. 

Alternately, each particle can be tracked individually in a Lagrangian frame of reference. In this 
formulation, position, velocity, mass, and temperature of each particle are updated at each time step based 
on mass, momentum and energy interactions with the gas phase and any surface reactions that may occur. 
This approach is typically used for very dilute particle concentrations with relatively small computational 
domains where the computational load is not unreasonable. This method was adopted by Egolfopoulos 
and Campbell (1999), and Zhou et al. (2002) in order to capture physical phenomena that could not be 
easily examined using the more traditional Eulerian approach. 

For this study, the Lagrangian formulation was adopted as the mass loading of particles was expected 
to be relatively low, and the consideration of larger particles made the continuum assumption required for 
the Eulerian approach appear questionable. Since these calculations are Quasi- ID, the calculations are 
additionally simplified by the lack of radial variation. This allows a single particle to be used 
computationally to represent all the particles at a given axial position, significantly reducing the number 
of particles that need to be analyzed at each time step. Particles are therefore inserted into the 
computational domain at the inflow boundary as sets of particles, with the number of particles in a set 
determined by the particle size and mass loading, for which only a single particle calculation must be 
performed per time step. The gas-particle interaction terms in the governing equations are then the results 
of the single particle calculation for a particle set at a given axial location multiplied by the number of 
particles in that particular set. This is particularly helpful for very small particles, where the number of 
particles can become very large even at small mass loadings. Once the gas-particle interaction terms have 
been calculated for a given set of particles, these values are linearly interpolated onto the gas-phase fixed 
grid for used in the gas-phase solution process. 

After a particle (or set of particles, hereafter referred to as though it were a single particle) is injected 
into the gas-phase flow with an initial mass, position, velocity, and temperature, a series of relations for 
spherical particles are used to determine the particle drag, heat transfer, and surface reaction rate that is to 
be applied to the particle. 

The first gas-particle interaction is aerodynamic drag. For laminar flow, White (1991) provides the 
following drag relation for a sphere. 


Fp 2 *~DPg u g 



(30) 


Many empirical relations have been proposed for the drag coefficient of a sphere. Carlson and 
Hoglund ( 1 964) and Crowe ( 1 967) developed correlations, based on the data available at the time, which 
have been used frequently for particle flow simulations. More current efforts, such as that of Igra and 
Takayama (1993) have provided more accurate expressions over a specific flow properties range based on 
detailed shock tube data. The drag coefficient used for this study is taken from Henderson (1976). This 
expression covers a very broad range of Reynolds number and Mach number with reasonable accuracy. 

As the relative velocity between the flow and the particles in this study could range from supersonic to 
near zero, and the static pressure from 80 atmospheres down to 0.1 atmospheres or less, it is important to 
cover the broadest possible range of the relative Mach and Reynolds numbers. The Henderson drag 
coefficient expression was developed from the extensive experimental data base compiled by Bailey and 
Hiatt (1971) through a series of ballistic range experiments and is accurate to within 16 percent across the 
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entire range considered. The expression is divided into three elements to cover the entire relative Mach 
number range. For a subsonic relative Mach number, 
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where the molecular speed ratio is defined as 


MSR=M reh l-L 


For a supersonic relative Mach number greater than 1 .75, 

l 


0.9 + 


0.34 

M re i 


+ 1.86 


M 


rel 


V Re ; -e/ J 


C D = 


„ 2 1.058 

2 + - + 


( rr, \ 


MSR 2 MSR 


K T sJ 


MSR 4 


f 


1 + 1.86 


M rel ‘ 

V ^ e re/ y 


(32) 


(33) 


In both the subsonic and supersonic equations, the relative Mach and Reynolds numbers are evaluated 
at free stream conditions. If the relative Mach number lies between 1.0 and 1.75, then the drag coefficient 
is taken as a linear interpolation between the value calculated using the subsonic Equation (31) with a 
Mach number of 1.0 and the value calculated using the supersonic Equation (33) with a Mach number of 
1.75, expressed as, 


C D {M } .ei , Re,- e/ ) = Qi (l .0, Re re/ ) + j (AT re/ - 1 .0\C D (l .75, Re,. e/ ) - Cp (l .0, Re,. e/ )) (34) 

The relative Reynolds number is evaluated at the actual free stream conditions for use in Equation (34). 
The gas-phase properties used in these equations are evaluated at the particle location, linearly 
interpolated between the nearest adjacent grid points. 

Once the force on the particle is calculated, the particle position and velocity can be updated by the 
simple relations, 
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Performing the energy balance on a large reacting particle in a convective flow has a number of 
complexities that need to be considered. Since the particle can be relatively large, a uniform particle 
temperature cannot be assumed, and an internal grid must be used to capture the temperature gradient 
within the particle. Since significant velocity slip is expected between the particle and the gas, convection 
must be modeled. Radiation losses must also be modeled since the gas surrounding the very hot, high 
emissivity particles is expected to be optically thin, and the test facility walls are water cooled. And, of 
course, since the particles are reacting, a large amount of heat is input to this system via combustion. 

If a thin shell at the outer surface of the particle is considered, as shown in Figure 3, an energy 
conservation Equation (37) can be written for that shell that encompasses each of the heat flows described 
above, and which allows for the accumulation or depletion of heat energy in the outer shell. 


Qaccum ~ Qcomb 9rad tfconv Qcond 


(37) 


Heat energy will then be conducted into or out of the particle interior over time at a rate controlled by 
the thermal diffusivity of the particle material. For the conduction, convection, and radiation terms, 
simple textbook relations are used, such as those found in Incropera and DeWitt (2002). 


Qcond ~ 


dT 
dr , 


(38) 


Outer Shell 



Figure 3. — Particle heat transfer and reaction model. 
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As with the drag coefficient, a number of empirical expressions for the Nusselt number have been 
developed for spheres over a variety of flow conditions, characterized by the Mach number and the 
Reynolds number. Once again, an expression that covers a wide range of these parameters is necessary for 
this study, including the case of supersonic relative velocity. In the case of supersonic relative velocity, 
the effect of the shock wave in front of the particle needs to be considered. Drake and Backer (1952) 
determined that it is best to correlate the Nusselt number with the conditions immediately behind the 
shock wave in line with the direction of travel. Since this is essentially a standing normal shock, simple 
textbook relations, such as those found in Saad (1985), can be used to calculate the pressure, temperature, 
and density behind the shock wave. After adjusting the thermodynamic and transport properties to the 
post-shock conditions, the post-shock Mach and Reynolds numbers can be calculated, and an appropriate 
correlation, such as that proposed by Fox, et al.(1977), given below, can be used. 
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For high speed flows, it is also necessary to calculate what is commonly referred to as the recovery 
temperature (White, 1988), given below for laminar flow. 
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(43) 


This value is then substituted into the convection heat transfer Equation (40), in place of the free stream 
gas temperature, to determine the convective heat transfer term. 

q C onv = ~h p 47 ir p 2 (T rec - T p ) (44) 


The recovery temperature is a modified form of the total temperature which approximates the temperature 
in the boundary layer surrounding the particle. This temperature is based on the free stream temperature 
and free stream relative Mach number, even in the case of a supersonic relative Mach number. 

The calculation of the combustion heat addition term, q com b, is not nearly so straight forward. A 
number of carbon particle combustion models have been developed, such as those of Blake (2002); Flurt, 
et al. (1998); Lee, et al. (1995); Makino and Law (1986); and Libby (1980). These models were typically 
developed to examine the oxidation of coal char. In addition to simplifications such as zero gas-particle 
velocity slip and uniform particle temperature, each of these studies was limited to combustion at lower 
pressures; typically atmospheric pressure. Flowever, for the problem under consideration in this study, a 
large range of pressures must be considered. Only two studies were found that looked at carbon 
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combustion at pressures considerably above atmospheric pressure in a convective flow, those of Moors 
(1998) and Hong (2000). The methodology of Hong was adopted for this study. 

Calculation of the combustion rate for solid carbon is somewhat complicated by the fact that the 
products of reaction shift as a function of the reaction temperature. As the carbon surface temperature 
increases, the reaction shifts from 


C(s) + 0 2 => C0 2 (45) 

to 

C(,)+^0 2 =>CO (46) 

Hong provides the following Arrhenius-type relation, taken from Hurt and Mitchell (1992), for the mole 
ratio of CO to C0 2 . 


-30000 

MR = moles co = 40000 g 1.9877;, (47) 

moles C o 2 

For convenience, this relation can be recast to give the ratio of moles of C0 2 produced by the surface 
reaction to the total moles produced. 


1 

M f = 

1 + MR 


(48) 


Two other convenient parameters can be defined at this point that will be useful for simplifying the 
particle combustion equations. These are the stoichiometric coefficient of oxygen consumed for each 
mole of carbon consumed, v 0 , and an intermediate variable, c. 

v 0 = 0.5(l + \|/) (49) 


% = 


\|/ + 1 


(50) 


Calculation of the combustion rate for a solid carbon particle is also complicated by the structure of 
the particle itself. For a typical carbon particle, the majority of the reactive surface area is made up of a 
labyrinth of micron to nanometer diameter pores internal to the particle. Oxygen must diffuse from the 
surface into the particle pore structure for significant reaction to take place. This process must be included 
in order to develop an accurate reaction rate model. This internal diffusion effect is captured through the 
use of an effectiveness factor, rj, defined by Hong as the ratio of the actual observed reaction rate 
including intra-particle diffusion effects to what the reaction rate would be if the oxygen concentration at 
the particle surface persisted throughout the interior pore structure of the particle. The effectiveness factor 
is calculated as follows. 


ft = /c 
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where 
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The Thiele Modulus, M T , is defined as 
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and the observed reaction order, m obs , is 
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The above expressions contain two kinetic rates for which Hong gives the following Arrhenius-type 
relations for non-porous graphite flakes, based on the experimental data of Ranish and Walker (1993). 

-51000 
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The constants contained in these kinetic rate expressions must be re-derived from experimental data for 
each individual type of graphite or coal char, as they are dependent upon the physical structure of the 
particles as well as their chemical composition. 

The Thiele Modulus expression also requires a value for the effective diffusion coefficient, D e . The 
effective diffusion coefficient is a function of the particle porosity, s (1.0 minus the ratio of the actual 
particle density to the density of a non-porous solid of the same material), the molecular diffusion 
coefficient of oxygen into the gas mixture at the local gas conditions, Dq,.v/ , and the Knudsen diffusion 

coefficient, D k (for diffusion in small pores where wall collisions dominate over inter-molecular 
collisions). Hong gives the following relation for the combination, 
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in units of cnr/s, as taken from Smith (1981). The value of Dq^m is available from the gas-phase 
portion of the Q1D code. 

In the formulation of Hong, a uniform particle temperature was assumed, so a single value of the 
effectiveness factor was used. In this study, it was necessary to use an approach that accounted for an 
internal particle temperature distribution in what is otherwise a global reaction model. The use of a mass 
averaged value of the effectiveness factor, equivalent to an internal surface area averaged value, was 
selected as an appropriate modeling approach. While no data was available to validate this modeling 
choice, this commonly used approach provides reasonable approximate values and reduces to the 
approach of Hong in the case of uniform particle temperature. The general behavior of the effectiveness 
factor as a function of temperature is shown in Figure 4, using the kinetic parameters for graphite flakes 
given above. 

In order to solve this complex system for the carbon burning rate, it is necessary to know both the 
overall particle oxygen consumption rate and the particle temperature profile. However, these values are 
inter-dependent, therefore requiring that the carbon burning rate be solved iteratively. Since the other heat 
flow terms from Equations (38), (39), and (44) depend on the particle surface temperature, and therefore 
the particle temperature profile, the energy balance given in Equation (37) must also be solved iteratively. 
This is done by running an inner iterative loop that corrects the particle surface partial pressure of oxygen 
(at the particle surface temperature supplied from an outer iterative loop) by comparing the 0 2 diffusion 
rate from the free stream and the total carbon consumption rate of the particle based on the graphite 
kinetics, taking into account the temperature dependent stoichiometric coefficients described above 
(Eqs. (45) to (47)). An outer iterative loop is run that corrects the particle surface temperature by 
comparing the outer shell carbon consumption rate calculated by an energy balance on the outer shell to 
the same rate calculated using the carbon consumption kinetic expression using the updated particle 
surface partial pressure of oxygen from the inner loop. 



0 500 1000 1500 2000 2500 


Temperature (K) 

Figure 4. — Effectiveness factor as a function of temperature. 
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The procedure for this solution is as follows 


1 . Assume a value for the particle surface temperature. This is typically selected as the value from 
the previous time step. 

2. Assume a value for the oxygen concentration at the particle surface. A value of 90 percent of the 
free stream oxygen concentration is used to start each iterative cycle. 

3. Calculate the global kinetic carbon consumption rate, r k , for the given values of the particle 
surface temperature, the particle interior temperature distribution, and the surface oxygen 
concentration. An intrinsic Langmuir expression is used to calculate the local volumetric 
consumption rate within each shell of the particle. This expression, adapted from Hong, is valid 
over a wide range pressures, temperatures and porosities. 


r k,i 



^int k\pjPo 2 ,sM c 
Sfot ^ L KpjPo 2 t S 


(59) 


The natural units for this expression is g/cm 3 -s, based on the particle volume, but it is 
typically given in terms of the particle external surface area, g/cnr-s. This transformation is easily 
done by multiplying the rate by the particle volume and then dividing by the external surface 
area. In either case, the units of k\ p and K p are mole/cm 3 -atm-s and 1/atm, respectively. The 
kinetic rates in Equation (59), k\ p and K p , are the same terms used in the Thiele Modulus 
expression. The interior, exterior, and total particle surface areas, S iat , S ext , and S tot , are used in the 
expression to proportion the amount of reaction that takes place on the particle exterior, and is 
therefore not affected by internal pore diffusion, versus that which occurs in the particle interior. 

It should be noted that the external surface area, S ext , is typically taken as the particle’s geometric 
external area multiplied by a roughness factor to account for surface irregularities. Hong suggests 
a value of 5 be used for the roughness factor. Using available data for porosity and pore diameter, 
and assuming a monodisperse pore size distribution within the particle, the particle interior 
surface area is readily calculated. 


16ner„ 
c P 

Oint — : 


1 pore 


(60) 


The local carbon consumption rate for each particle shell volume is calculated by multiplying r ki 
by the volume of the particle shell being evaluated. The global carbon consumption rate is then 
obtained by summing the individual shell rates. This value is then normalized by the particle’s 
geometric external surface area, so as to be comparable to the other rates used in the iterative 
procedure, which are expressed per unit external surface area. 

4. Next, the rate of carbon consumption based on diffusion of 0 2 from the free stream to the particle 
surface must be calculated. The mass transfer coefficient for 0 2 diffusing to the surface, recast in 
terms of carbon consumption, can be expressed as, 

M c D 0l „,Sh 

Kq — 

2 Kp R u TfUm V o 

where 

Tfilm = 0 - 5 ( t p +T g ) 


(61) 


(62) 
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and where Sh is the Sherwood number, which is calculated from Equation (42) by substituting the 
Schmidt number for Prandtl number. A Schmidt number of 0.74 is used for this study. Then the 
rate of carbon consumption based on oxygen diffusion to the surface is given by, 
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(63) 


or 


r ox = k d {Po 2 ,f s - p o 2 ,s) (^ - 0) (64) 

with the pressure terms in the above expressions given in atmospheres. 

Assuming that the reaction and diffusion processes are in a quasi-steady state, the carbon 
consumption rate from Equation (63) or (64) and the rate of carbon consumption calculated from 
the kinetic expressions above should be equal. If they are not initially equal, the surface partial 
pressure of 0 2 Po 2 ,s > is iterated upon to bring these two rates into equilibrium. 

5. The carbon consumption-based 0 2 diffusion and 0 2 consumption rates are now in equilibrium at 
the particle surface temperature assumed in step 1 . An energy balance of the particle shell using 
the terms shown in Figure 3 is then used to adjust the particle surface temperature as follows. 
First, (37) is rearranged to solve for q com b. 


({comb ~ q accum + qrad qconv + qcond 


(65) 


where 


q accum 


47tp iS 


3 


A r 


C S (T P M -TJ). 


( 66 ) 


Then we can write a carbon consumption rate expression for the outer shell based on the rate that 
would be required to produce q comb . 


„ qcomb 

e- (l- V )A// C o+vAtfco 2 


(67) 


If r e is not equal to the kinetic rate for the outer shell calculated using Equation (59) within a 
given convergence criteria, T p is adjusted upward or downward, as appropriate, and the procedure 
is repeated starting at step 2 until convergence is achieved. 


Now that the carbon consumption rate and particle outer shell energy balance have been calculated, 
the next step is to update the particle mass to account for the mass lost due to the interior and exterior 
surface reactions. Mass can be lost from both the exterior surface of the particle, reducing the particle 
diameter, and from the pore surfaces internal to the particle, increasing the diameter of the pores. To 
simplify the analysis for this study, it was decided to adjust the particle mass by increasing the pore 
diameter within the particle. This approach was selected because of the high porosity of the graphite 
material being studied and the difficulty in establishing a proper mass decrease to diameter decrease ratio. 
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According to Hong, the relationship between mass decrease and particle diameter is generally determined 
experimentally, and this ratio changes during the combustion history of a particle and with ambient 
conditions. By changing the pore diameter only, a computational efficiency is obtained in that the particle 
grid does not have to be adjusted after each iteration. 

Once this has been done, the source terms for overall mass, momentum (including the momentum 
associated with the mass lost by the particle), energy, and species (CO, C0 2 , and 0 2 ) are calculated for 
use in Equations (1), (2), (3), and (6). These source terms are multiplied by the number of particles in 
each particle set, and then linearly interpolated onto the two adjacent grid points for use in the gas phase 
governing equations. 

Lastly, the internal particle temperature profile is updated using the transient, 1 -D conduction 
equation in radial coordinates. 


j__a/ 

r 2 dr \ 


K s r 
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Ps C s St 


( 68 ) 


The simple, 2 nd -order accurate Crank-Nicolson semi-implicit method (Tannehill, et al., 1997) is used 
to solve the conduction equation. This method was selected because it is unconditionally stable and does 
not limit the time step used in the overall calculation. 


2.7 Gas Phase Numerical Solution Method 

The gas-phase governing equations (Eqs. (1), (2), (3), and (6)) are cast in the following form for 
solution 
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Two different schemes are available in the code for solving this set of equations. The first scheme is 
the standard MacCormack predictor-corrector scheme, as described in Anderson (1995). This solution 
method is 2 nd order accurate in space and time, and follows the following steps. 


Predictor Step 


1 . The values in the F and J arrays are first calculated from the values contained in the U array from 
the previous time step. 

2. Equation (69) is then used to obtain the time derivative of the U array using a forward difference 
in x. 

3. The U array is then updated to an interim value by the equation, 


U =UJ 


+ 


d dU N 

V St , 


A t + % 


(73) 


where j refers to the time step and x is a Total Variation Diminishing (TVD) artificial viscosity 
term that will be discussed later. 


Corrector Step 


1 . The values in the F and J arrays are updated from values contained in the U array. 

2. Equation (69) is again used to obtain the time derivative of the U vector using a backward 
difference in x. 

3. An average of the two time derivatives is calculated as follows, 
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(74) 


4. Lastly, the U array is updated to the next time step by the following expression, 


UJ +1 =UJ +\ 


dU 
dt , 


A t + X 


(75) 


avg 


This scheme can be easily improved to 4 th order accuracy in space by making the following 
substitutions (Georgiadis 2000), with k representing the spatial grid coordinate. 


1 . Calculate the time derivative of U from predictor step 2 by the expression, 

^T = -Tr(- 1F k + 8F* + 1 - F k+2 ) + J (76) 

dt 6Ax 

2. Calculate the time derivative of U from corrector step 2 by the expression, 

^r = ~TT [ 1F k - 8 Fk-I + F k -1 ) + j ( 77 ) 

dt 6 Ax 

The user can select either scheme within the code. 
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Numerical viscosity, /, which is generally required when using a MacComiack type scheme 
(Anderson 1995), is combined with a TVD routine in an easily implemented method developed by Davis 
(1987). The basic form of the TVD artificial viscosity term is 


K+ k+\ 
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The leading variable in this expression is defined as 
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where the Courant number, v„ is defined by 
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and the eigenvalue at each point is defined as 
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The expression for the flux limiter ©(r*) is given by 


eH= 
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(83) 


where r + and r are defined by the following inner products, 


AU k-V’ AU k+v 


w k+ y 2 ’ AU k+ y 2 


(84) 
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AU k-y 2 ’ AU k+ y 


AU k-y 2 ’ AU k-y 


(85) 


and AU represents a backward difference in ,x. 

The gas-phase solution method also requires boundary conditions. These boundary condition may 
need to represent subsonic, choked, or supersonic flows, depending on the simulation. For supersonic and 
choked inflow conditions, the velocity or Mach number is set along with either the static or total pressure 
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and temperature. For supersonic and choked outflow conditions, simple extrapolated boundary conditions 
are used. Subsonic boundary conditions are inherently more challenging, since information can travel 
both upstream and downstream. The methodology of characteristic boundary conditions of Poinsot and 
Lele (1992), valid for both the Euler equations as well as the Navier-Stokes equations, has been adopted 
for these conditions. The boundary condition can be pre-defined, or can be selected by the code as a 
function of the flow conditions as the calculation progresses. 

The computational grid used for these calculations is generally uniform in the axial coordinate x, with 
grid spacings ranging from approximately 4 mm down to 1 mm. Since the calculations are quasi- ID, the 
number of grid points required for any given calculation is typically reasonable, from 300 up to 4000. 
Departures from a uniform grid are pointed out where they occur. 

The time step used in these calculations was allowed to vary during the course of the computation to 
the largest value that did not exceed any of three limits put in place to ensure stability and accuracy of the 
solution. The first limit is the standard Courant, Friedrichs, and Lewy (CFL) stability limit for explicit 
numerical methods (Tannehill, et al., 1997). A time step limit of 50 percent of the maximum time step 
allowed by the CFL criterion was used to ensure stability. A gas-phase chemical reaction rate limit to the 
time step was also found to be necessary to prevent a chemical heat release “runaway” under some 
conditions. A limit of 1 percent change to any species mass fraction at any grid point was typically used. 
Lastly, it was found during validation testing of the code that a time step limit based on the rate of change 
of the outer shell temperature of the particles was necessary to maintain accuracy. This limit was typically 
set to 0.02 K temperature change in the outer shell of any particle. The smallest of these three limits is 
selected at each time step as the value to be used for the succeeding time step. 

3.0 Validation 

For the purposes of validating the above described Q1D 2-phase code, no single representative test 
case of known solution was found. Instead, different aspects of the code were validated individually using 
available solutions. Individual physical process validations were performed in addition to the higher-level 
validations included in these sections. These physical process validations included chemical equilibrium 
calculations, shock tube calculations, individual particle fixed outer temperature boundary condition 
calculations, individual particle constant heat flux calculations, rocket performance calculations, steady 
state nozzle flow calculations, and various carbon combustion calculations. Additionally, many test cases 
were run to determine maximum acceptable time steps and grid sizes and to determine appropriate 
convergence criteria for a wide range of problems. 

3.1 Gas Phase Flow and Chemical Kinetics Validation 

As a first step, the gas-phase reacting flow capability was validated using a 1-D detonation 
simulation. The full kinetic mechanism was exercised by using a reactant mixture of 10 percent CO, 

3 percent H 2 , 32 percent 0 2 and 55 percent N 2 by mass. As shown in Figure 5, the detonation was initiated 
by a small region of high pressure, high temperature FL at the closed end of a tube filled with reactants at 
300 K and 100 kPa. The ignition region is initially at 3000 K and 4 MPa. The detonation wave stabilized 
after traversing the first 7 mm of the computational domain, so the calculation was terminated at that 
point. A time step of approximately 1.5x10 10 s was used. 

These results were then compared to a calculation obtained from the 1 -dimensional ZND code 
developed by Shepherd (1986). To run the ZND code, the desired kinetic mechanism is first input into the 
CHEMK1N code, described by Kee et al. (1993), that provides the kinetic and thermodynamic properties 
required. Next, the Gordon and McBride (1994) CEA code, or any other equilibrium detonation code, is 
used to determine the detonation (Chapman-Jouguet) velocity. With these inputs, ZND calculates the 1-D 
detonation structure by first assuming a normal shock based on the detonation velocity and quiescent 
reactant conditions and then reacts the mixture via the kinetic mechanism and gas properties from 
CFIEMK1N using the post-shock conditions as the reaction initial conditions. ZND performs these 
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calculations in the reference frame fixed to the detonation wave, shown in Figure 6, so the velocity must 
be transformed back to the stationary reference frame by adding the detonation wave velocity. All of the 
static properties are unaffected by the change in reference frames. 

Figure 7 shows the convergence of the detonation wave speed upon the equilibrium value calculated 
using the CEA code. The wave speed at the termination of the calculation is approximately 1 percent 
below the CEA value, which was judged to be sufficiently converged for the purposes of this validation. 



Figure 5. — Detonation case geometry. 
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Figure 6. — Detonation calculation reference frames. 
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Figure 7. — Detonation wave speed convergence. 


Figures 8 to 15 show the values of pressure, temperature, velocity, and H 2 0, OH, C0 2 , NO, and N0 2 
mass fractions calculated by Q ID at the end of the computation against those calculated by ZND using 
the same kinetic mechanism given in Appendix B. Since the results from ZND begin after the leading 
shock wave of the detonation front and therefore do not include part of the induction zone, leading to a 
difference in total induction time between the two codes, the two sets of results were aligned by matching 
the midpoints of the initial temperature rise region of the temperature profiles shown in Figure 9. Once 
the results were thus synchronized, the ZND and Q1D results generally agree quite well. The largest 
differences are observed in the region of the shock, particularly noticeable in the pressure and velocity 
profiles. This is caused by the chemical reactions and consequent heat release initiating within the shock 
wave in the Q1D model, whereas in the ZND model, it is assumed that the chemical induction time starts 
after the passing of the shock wave, as mentioned above. This difference in the time at which reactions 
start is evident in the H 2 0, OH, and C0 2 profiles. The nitrogen reactions are somewhat slower, and 
therefore do not reflect this difference as significantly. The quicker rise in NO concentration shown in 
Figure 14 for the Q1D code is most likely due to differences in the thermodynamic and transport property 
data sets used by the two codes. Overall, the good agreement between these two results for the detonation 
case indicates that the gas phase fluid dynamic and kinetic models and solvers within the Q1D code are 
working properly. 
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Figure 8. — Comparison of pressure profiles between Q1D and ZND codes. 
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Figure 9. — Comparison of temperature profiles between Q1D and ZND codes. 
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Figure 10. — Comparison of velocity profiles between Q1D and ZND codes. 
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Figure 1 1 . — Comparison of H 2 O mole fraction profiles between Q1 D and ZND codes. 
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Figure 12. — Comparison of detonation OH mole fraction profiles between Q1D and 
ZND codes. 
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Figure 13. — Comparison of detonation CO 2 mole fraction profiles between Q1D 
and ZND codes. 
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Figure 14. — Comparison of detonation NO mole fraction profiles between Q1D and 
ZND codes. 
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Figure 15. — Comparison of detonation NO 2 mole fraction profiles between Q1D 
and ZND codes. 
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3.2 Particle Drag and Heat Transfer Models and Gas-Solid Momentum and Heat 
Exchange Models Validation 

The next step in validating the Q1D code is to verify that the particle momentum exchange and 
convective heat transfer models are working properly. This can be accomplished by studying a gas flow 
interaction with chemically inert particles. Maxwell and Seasholtz (1974) examined the behavior of 
micron-size inert particles passing through a Mach 1.6 normal shock wave for the purposes of 
determining the tracking behavior of these particles used in Laser Doppler Velocimetry (LDV). The 
physical situation modeled in the study was that of a constant area duct with a supersonic, 289 K static 
temperature, 1 atmosphere static pressure, particle laden inlet air flow and a standing normal shock wave 
interior to the duct. The particles were assumed to be at the gas static temperature and velocity up to the 
shock wave, and to occupy negligible volume. For the analysis, the post-shock conditions, as calculated 
using the standard Hugoniot relations, were used as the initial conditions for the gas, while the pre-shock 
temperature and velocity were used as the initial conditions for the particles. Since only small particles, 
less than 5 pm diameter, were considered for this study, uniform particle temperature was assumed 
throughout. Variations in particle diameter, particle density, gas and particle inflow velocity, particle 
specific heat and particle/gas mass loading ratio were examined. Constant thermodynamic and transport 
properties were assumed for the air and the particles, but the exact values were not given. Because of this, 
the calculations had to be replicated using the following equations adapted from Maxwell and Seasholtz 
(1974), as the values assigned for the thermodynamic and transport properties can significantly affect the 
results. 

A force balance on an individual particle yields the following relation, 
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where the absolute value of the gas-particle velocity difference is used as shown to maintain the correct 
sign for the particle velocity spatial derivative. Likewise, an energy balance on an individual particle 
yields the following equation, 
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Equating the momentum and energy transfers between the two phases, with a defined as the ratio of the 
particle to gas mass flows, yields the following two conservation equations. 
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In order to solve this system of differential equations, it is also necessary to use the following relation for 
the rate of change of pressure, generated by taking the spatial derivative of the ideal gas law and 
transforming the density term into a velocity term through the expression for continuity in a constant area 
duct. 


8P _ P ( dT g T g dV g ' 
dx T g ^ dx V g dx ^ 


(90) 


The particle drag coefficient and Nusselt number correlations used in the above equations are taken from 
Maxwell and Seasholtz, and differ only slightly in value from the expressions used in the Q1D code in 
this flow regime. These equations were then solved using a simple explicit space-marching routine, with a 
step size of 1 x 1 0 6 m. 

In order to fully validate the two-phase capability of the Q1D code, a particle/gas mass ratio, a, of 0.2 
was selected for comparison with a solution using the analysis method of Maxwell and Seasholtz. At this 
mass ratio, both the particles and the gas respond significantly to the presence of the other phase. A 
particle diameter of 2 pm and a particle/gas specific heat ratio of 1.125 was used in both analyses. The 
Q1D simulation was run by first establishing a steady, subsonic gas flow in the constant area duct that 
approximated the post-shock conditions, as determined from the standard Hugoniot relations. A total 
pressure, total temperature subsonic characteristic inflow boundary condition and a fixed static pressure 
subsonic characteristic outflow boundary condition, both based on the methodology of Poinsot and Lele 
(1992), were used to establish this initial solution. The particles were then introduced into the flow at the 
desired mass fraction from a fixed position in the computational domain as the solution continued to run. 
Each set of particles was inserted at the pre -shock static temperature and flow velocity. The solution was 
then allowed to run out to a new steady state solution. As the gas-solid interactions caused the gas mass 
flow to shift from the initial result, the particle mass flow rate was continuously adjusted to maintain a 
particle/gas mass ratio of 0.2. The gas phase inlet velocity, static temperature, and static pressure also 
shifted from the initial gas-phase solution due to the interaction with the particles. These new static/total 
ratios were then characteristic of a slightly different shock Mach number. Therefore, the analytical 
solution was run at the gas flow inlet conditions of the converged numerical solution for comparison. 

Figures 16 and 17 show excellent agreement between the Q1D simulation and the M ax we I l/S cas h o I tz 
analysis gas temperature and gas velocity from the point of particle introduction to the end of the 
computational domain. Figures 18 and 19 show the corresponding results for the particle “stream”, also 
with excellent agreement between the two methods. 

These results are consistent with the expected behavior of the two phases as given by Maxwell and 
Seasholtz. The slight differences between the results of the two methods observed in Figures 1 6 to 19 are 
due to a combination of the different particle drag coefficient and Nusselt number correlations used in the 
different calculations, and the use of varying thermodynamic and transport properties in the Q1D code 
versus the fixed values used in the MaxwelFSeasholtz analysis. 
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Figure 16. — Gas temperature profile from Q1D calculation from point of particle 
insertion compared to gas temperature profile from Maxwell/Seasholtz 
analysis. 
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Figure 17. — Gas velocity profile from Q1D calculation from point of particle 
insertion compared to gas velocity profile from Maxwell/Seasholtz 
analysis. 
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Figure 18. — Particle temperature profile from Q1D calculation from point of 
particle insertion compared to particle temperature profile from 
Maxwell/Seasholtz analysis. 
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Figure 19. — Particle velocity profile from Q1D calculation from point of 
particle insertion compared to particle velocity profile from 
Maxwell/Seasholtz analysis. 
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3.3 Particle Internal Heat Transfer, Heat Release and Supersonic Drag/Convection 
Validation 

Next, the results of Baek (1985), summarized also in Sichel, et al. (1985), were used to validate the 
particle internal heat transfer/heat release modeling in the Q1D code, as well as the supersonic relative Mach 
number drag and convective heat transfer correlations. Baek performed a series of experiments examining 
the shock ignition of various types and sizes of particles, which he then also examined using a 
computational model and asymptotic analysis. The most comprehensive set of results reported by Baek were 
obtained using 53 pm coal particles subjected to a Mach 4.8 traveling shock wave. The particles were 
injected into a shock tube containing pure oxygen at 1/3 of an atmosphere pressure and 295 K in the driven 
section. The mass of the particles was very small compared to the mass of gas present, so the particles were 
assumed to have no effect on the gas flow. The particles were also assumed to be of uniform size and 
spherical. The coal particles were characterized as having the following physical properties. 


Density: 

Thermal Conductivity: 
Heat Capacity: 

Heat of Combustion: 
Internal Surface Area: 


1200 Kg/m 3 
0.887 J/m-s-K 
987.4 J/Kg-K 
35812 KJ/Kg 
426000 m 2 /Kg 


Under the conditions tested, the particle consumption rate, in g/cnr-s (based on the particle external 
surface area), was assumed to follow a simple Arrhenius expression of the same form as shown in 
Equation (91), with an activation energy of 35,700 cal/g-mole and a pre-exponential factor of 
87,100 Kg/nrt-s-atm. Baek’s numerical single-particle combustion model was found to match the 
experimental data well over a range of shock Mach numbers and particle sizes. 

In the Q1D code, this scenario was modeled by looking at a single 53 pm coal particle, with the 
physical properties listed above, injected at zero velocity into a steady flow at the post-shock conditions 
as calculated from a standard ideal shock tube calculation, such as that found in Saad (1985). A particle 
porosity of 0.32 with 2 nm diameter pores was found to give the required internal surface area. The two 
Arrhenius rate expressions required for the Hong particle combustion kinetic rate, Equation (59) (restated 
below), were defined as shown in Equations (9 1 ) and (92) to match the behavior of the simple kinetic 
expression of Baek. The Q1D code was then run until ignition occurred. 


= 
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-35700 

k\ p - 3.09xl0 9 e r,,t p j (91) 

-30300 

K p =13.4e R » T pJ (92) 

The results from the Q1D code are compared to the results from Baek (1985) in Figures 20 to 22. The 
Baek results were recalculated using Baek’s original code, a listing of which was included in the report. 
Figure 20 shows the gas-particle relative Mach number from just after particle insertion up until the time 
at which the particle ignited in the Baek calculation. The small difference in relative Mach number 
reflects differences in the thermodynamic properties between the two codes, resulting in a slightly 
different sound speed in the gas. Other than this slight offset, the Mach number decay is nearly identical, 
which is to be expected, as Baek also used the Henderson (1976) particle drag correlations, given in 
Equations (31) to (34). 


N ASA/TM— 20 10-216765 


34 



Temperature (K) 







Figure 21 shows a time series of particle internal temperature profiles as calculated using Baek’s 
code. The particle surface temperature is at the left hand side of the plot, with the temperature falling off 
with distance into the particle. The lowest temperature profile corresponds to 2 psec after the passage of 
the shock wave, and each succeeding temperature profile represents the passage of another 2 psec of time. 
Ignition occurs just prior to the 14 psec (uppermost) profile. 

Figure 22 shows a similar plot of results from the Q1D code. The overall ignition process took 
approximately 4 psec longer in this simulation, due to a combination of small differences in the 
thermodynamic properties, transport properties, and heat transfer relations used in the different codes. 
While this approximately 25 percent increase in ignition time seems significant, it actually represents a 
very small difference in particle heating rate. The continual decrease in the convective heating rate of the 
particle due to the decrease in the relative velocity between the gas and the particle and the decreasing 
temperature difference leads to the particle surface approaching the ignition temperature asymptotically. 
This behavior translates to a high sensitivity of the ignition time to the heating rate, since it could take a 
significant portion of the ignition time to heat the last few percent of the required temperature rise. 

Even with this longer ignition time, the Q1D result is within the scatter of the experimental data for 
ignition delay given by Baek for this condition. Taking this into consideration, the Q1D code results 
compare quite favorably with the analysis and experimental results of Baek. 



Figure 22. — Internal temperature profiles at 2 psec intervals of particle heating up to 
ignition from Q1D code. 
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4.0 Generation of Graphite Kinetic Constants 

The graphite combustion rate expression. Equation (59), repeated below, requires that the constants in 
the kinetic rate expressions, Equations (55) and (56), also repeated below, be derived from experimental 
data for the given material and particle shape. 
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The above expressions contain the two sets of kinetic rate constants given by Elong (2000) for non- 
porous graphite flakes, based on the experimental data of Ranish and Walker (1993). While the data used 
to generate this data covers a broad pressure range (from 0.1 MPa up to 6.4 MPa) over a solely kinetically 
controlled temperature range (733 to 842 K), the graphite flakes used were significantly different in 
composition, structure and shape than the graphite spheres of interest in this study. The Ranish and 
Walker graphite flakes averaged 30 pm in diameter and 0.5 pm thickness. They were also a high purity, 
non-porous natural crystalline graphite. In the HTF facility described earlier, the graphite used in the 
facility heater and the piping liners is Union Carbide PGX graphite (GrafTech, 2005). PGX is an artificial 
molded graphite material formed using lampblack particles as large as 0.8 nun held together with a coal- 
tar pitch based binder. These materials are mixed together, poured into a mold and graphitized by baking. 
The resultant material has a density of 1730 Kg/m 3 (Graftech, 2005), which translates to a porosity of 
0.2345 (based on a non-porous graphite density of 2260 Kg/nr 3 ) with an average pore diameter of 6 pm 
(Pierson, 1993). Observations made during test facility operations indicate that roughly spherical particles 
are to be expected. Because of these differences, an alternate data set from which to generate the kinetic 
constants was needed. However, no available data set was found in the literature for an appropriate 
material that covered the pressure and temperature range of interest. In particular, no data sets including 
data for pressures above 1 atmosphere were found. Because of this, the Hong kinetic constants were used 
as a starting point, since they included the effect of pressure, and then the constants in the k Xp expression 
were modified to fit the available ambient pressure data. 

The frequently referenced graphite sphere combustion rate data of Tu, et al. (1934) was one possible 
data set, but upon examination, it was found that the graphite used was not sufficiently characterized so as 
to provide all the inputs needed for the combustion model of Hong. A number of other graphite oxidation 
data sets, such as that of Okada and Ikegawa (1953), were obtained for cylindrical graphite electrodes, but 
not for spheres. The only applicable data was that taken by Froberg (1967) in a set of experiments 
modeled after the work of Tu, et al. (1934), but for which significantly more graphite characterization 
information was available. The graphite used by Froberg had a density of 1650 Kg/nr 3 and an average 
pore diameter of 5 pm, both similar to the PGX graphite used at HTF. In the set of Froberg’ s experiments 
that were used to generate the graphite rate constants for this study, large graphite spheres (-12.6 mm 
diameter) were exposed to heated air streams with air temperatures ranging from 8 1 3 to 1 1 96 K. The air 
heater inflow velocity was held constant at 1.0 cm/s (at standard temperature and pressure), which 
resulted in an increasing air velocity as the air temperature was raised. The graphite sphere was allowed to 
come to equilibrium before combustion rate data was taken, which resulted in graphite surface 
temperatures as high as 1311 K for the highest temperature air flow. The combustion rate was directly 
measured via sphere weight loss. 
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For the Q1D simulation of this configuration, thermal conductivity and heat capacity as a function of 
temperature were also needed for the graphite. As this data was not available from Froberg for the 
specific graphite studied, estimates were made based on literature information. The thermal conductivity 
of PGX is 120 W/m-K at room temperature (Graftech, 2005) and decays exponentially as temperature 
increases. It was assumed that the Froberg graphite had the same room temperature thermal conductivity 
as the PGX. The thermal conductivity decay rate, as well as graphite specific heat as a function of 
temperature, are given graphically by Pierson (1993), and were used to generate 3 ld order polynomial 
curve fits of these properties for use in the Q1D code, as shown in Figures 23 and 24. For those cases 
where an internal particle temperature profile was used, mass-averaged values for heat capacity and 
thermal conductivity were used. 




Temperature (K) 

Figure 24. — PGX graphite heat capacity. 
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Figure 25. — Q1 D kinetic constant fit to Froberg data. 



With this information, Q1D simulations of a number of data points from Froberg were made, and the 
kinetic constants adjusted until the simulations and the experimental results were in reasonable 
agreement. The resultant expression for k Xp is given in Equation (93), which represents a factor of 1 1.3 
increase in the pre -exponential factor and a 24 percent increase in the activation energy over the 
expression from Hong (2000) for non-porous graphite flakes given in Equation (55). 

-63100 

k\ p = 7. 1 1 x 10 9 e RuT pj (93) 

A graphical comparison between the Froberg data and the adjusted kinetic model is shown in Figure 25, 
based on the sphere surface temperature. Using this adjusted expression for k\ p , PGX graphite can be 
simulated by adjusting the particle density, pore diameter, and particle size to the appropriate values. 


5.0 Model Facility and Scramjet Simulation Geometry and Flow Conditions 

5.1 Test Facility Components 

In the introduction, two different hypersonic propulsion test facility storage-heater types were 
described. The first utilized ceramic pebbles heated by flowing combustion exhaust gases through the 
bed. The second consisted of inductively heated drilled graphite blocks. The first type is found in the 
GASL Leg IV facility (Roffe, et al., 1997), while the second type is used at the NASA Plumbrook HTF 
facility (Perkins, et al., 1998). Schematics of these facilities were presented in Figure 1 and Figure 2. 
Rather than attempting to simulate both of these existing facilities for this code demonstration, a single 
facility geometry will be used to allow for a direct comparison of results for non-reacting and reacting 
particles. As the goal of this study is not to determine the state of particulates entering the combustor of a 
particular facility at a specific test condition, but instead to parametrically determine the range of 
particulate states that could be encountered, and the dominant physical processes that determine these 
states, the use of a generic facility geometry does not detract from the objectives of the study. 
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300 cm 300 cm 



Figure 26. — Scaled axisymmetric geometry for facility settling chamber/mixer and nozzle. 

The selected generic facility geometry, shown in Figure 26, consists of a 3 m long low speed 
(~ Mach 0.08 gas flow velocity) mixing section followed by a 3 m long Mach 5 exit freejet nozzle. The 
scale of the mixer and nozzle sections is slightly less than that of the NASA HTF facility, and 
considerably larger than that of GASL Leg IV facility. The nozzle contour was generated using a spline 
fit to the inlet, throat, and exit axial locations and diameters, estimated from the facility descriptions given 
in the above cited facility references. The same nozzle geometry is used for testing at simulated flight 
Mach numbers of 5, 6, and 7 by adjusting the mixer section inflow total pressure and total temperature, 
and by adjusting the flow deflection angle and length of a pre -compression plate that will be described in 
detail in the next section. The inflow conditions for the Mach 5, 6 and 7 simulated flight conditions are 
given in Table 1, along with the representative flight dynamic pressures used to calculate the facility 
conditions, and include corrections for real gas effects taken from NASA compressible flow tables (Ames 
Research Staff, 1953). 


TABLE 1.— MIXER SECTION INFLOW TOTAL PRESSURE. TOTAL TEMPERATURE AND 
DYNAMIC PRESSURE FOR MACH 5, 6 AND 7 SIMULATED FLIGHT CONDITIONS 



P total 

(KPa) 

T total 

(K) 

Dynamic pressure 
(KPa) 

Mach 5 

3082 

1219 

89 

Mach 6 

4498 

1652 

54 

Mach 7 

8247 

2174 

47 


5.2 Scramjet Inlet and Isolator Sections 

After the facility components are defined, it is necessary to define the geometry of the scramjet model 
to be simulated. In practice, the entire scramjet inlet is not used in ground testing due to its length. The 
multiple external compression ramps are instead replaced by a single pre-compression plate that provides 
the same inflow condition at the engine inflow plane. Flowever, it is still necessary to define the full inlet 
geometry in order to determine the appropriate pre-compression plate geometry. Figure 27 shows a 
generic, unsealed scramjet inlet and isolator, along with a short section of the combustor. Table 2 defines 
the geometric parameters shown in the figure for a Mach 7 design point inlet/isolator system. The values 
upstream of the engine throat were determined by first selecting the throat height, H l5 to set an overall 
engine scale consistent with the generic facility size, and then iterating on the compression ramp angles 
and lengths such that all of the oblique shocks converge on the cowl lip, and the cowl lip oblique shock in 
turn impinges on the shoulder that defines the inlet throat, as shown in Figure 27. This design was further 
constrained by the typical design practices of keeping the compression ramp turning angles approximately 
equal to minimize overall total pressure loss, and of targeting a throat Mach number of approximately 
one -half of the free stream Mach number. The isolator section was sized using the methodology described 
by Heiser and Pratt (1994) for constant area isolators, with an isolator pressure ratio of 2. The slight 
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divergence angle within the isolator section accounts for boundary layer growth that would otherwise 
restrict the flow in the duct. A series of linked spreadsheets were used to automate this design process. 




L \ 




Figure 28. — Scaled inlet/isolator at Mach 7 design point. 


Figure 28 shows a scaled representation of the inlet/isolator and the shock system at the Mach 7 
design point based on the geometric parameters given in Table 2. Figures 29 and 30 show the shock 
structure moving out from the cowl lip as the flight Mach number is decreased to Mach 6 and Mach 5, 
respectively. Table 3 then provides the static pressure, static temperature, and Mach number at key points 
through the geometry for the selected Mach 5, 6 and 7 flight conditions. It should be noted that the 
conditions at the isolator entrance are a result of the combination of the cowl lip oblique shock with an 
expansion fan emanating from the inflection point at the inlet throat , which is not shown in Figures 27 
to 30. 


TABLE 2.— SUMMARY OF MODEL SCRAMJET INLET/ISOLATOR 
GEOMETRY AT MACH 7 DESIGN POINT 


L { = 353.2 cm 

H l = 6.0 cm 

O 

SO 

II 

c£l 

L 2 =120.2 cm 

P 

II 

00 

O 

02 = 8 . 2 ° 

L 3 = 99.6 cm 

P 

II 

O 

03 = 10.4° 

= 33.2 cm 

cx 3 = 26° 

04=13.2° 

Z 5 = 547.6 cm 

a 4 = 2 ° 
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Figure 29. — Scaled inlet/isolator with shock locations corresponding to Mach 6 operation. 



Figure 30. — Scaled inlet/isolator with shock locations corresponding to Mach 5 operation. 


TABLE 3.— SUMMARY OF FLOW FIELD CONDITIONS FOR MODEL SCRAMJET 



Mach 5 

Mach 6 

Mach 7 | 


P 

1 static 

T 

J static 

Mach 

P 

± static 

T 

J static 

Mach 

P 

± static 

T 

J static 

Mach 


(Pa) 

(K) 

no. 

(Pa) 

(K) 

no. 

(Pa) 

(K) 

no. 

Free stream 

5024 

215 

5 

2110 

227 

6 

1370 

230 

7 

Along 1 st ramp 

12600 

285 

4.2 

6200 

319 

5.0 

4670 

342 

5.6 

Along 2 nd ramp 

30100 

372 

3.6 

16900 

434 

4.1 

14200 

483 

4.6 

Along 3 rd ramp 

63500 

463 

3.0 

39300 

555 

3.5 

35800 

634 

3.8 

Isolator inlet 

60400 

512 

2.8 

37400 

636 

3.1 

34300 

751 

3.4 

Isolator exit 

120800 

798 

1.8 

74800 

1000 

2.0 

68600 

1237 

2.2 


Having established the scramjet inlet/isolator geometry and the flow conditions at each step through 
those components, the test configuration to be used can be defined. To shorten the scramjet geometry, the 
first two external compression ramps from Figure 27 are replaced by a single pre-compression plate, 
leaving the last external compression ramp, the cowl, and the isolator unchanged. In order to create the 
same conditions along the last external compression ramp before the throat as those achieved in the 
original configuration from Figure 27, the geometry has to be rotated slightly counter-clockwise, as 
shown in Figure 3 1(a). This rotation compensates for the change in the flow direction coming at the last 
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external compression ramp in the shortened configuration relative to that in the original configuration. 
Each of the geometric parameters shown in Figure 3 1(a) needs to be changed for each test Mach number 
to be able to match the flow conditions approaching the last external compression ramp calculated for the 
original configuration. The required values as a function of flight simulation Mach number are given in 
Table 4. The flow conditions remain the same along the last external compression ramp and within the 
isolator, but a new set of flow conditions are created along the precompression plate, which are given in 
Table 5. Figure 31 (b) shows the test article from Figure 31(a) schematically “installed” in the model test 
facility, indicating the full computational domain. 


Final External 

Compression Ramp Isolator 



(a) 

O 2 & N 2 Injection 



Figure 31 . — (a) Truncated configuration for wind tunnel testing with forebody simulation (precompression) 
plate and (b) truncated configuration with forebody simulation plate shown schematically with facility mixer 
and nozzle indicating computational domain. 
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TABLE 4.— TEST ARTICLE CONFIGURATION FOR MACH 5, 6 AND 7 TEST CONDITIONS 



Mach 5 

Mach 6 

Mach 7 

p 5 

9.3° 

9.3° 

10.1° 

a 5 

15.0° 

9.5° 

4.5° 

a 6 

2.0° 

7.5° 

12.5° 

U 

100.6 cm 

99.1 cm 

84.0 cm 


TABLE 5.— FLOW CONDITIONS ALONG THE PRE-COMPRESSION PLATE 
FOR MACH 5, 6, AND 7 TEST CONDITIONS 



Mach 5 

Mach 6 

Mach 7 

P static 

23700 

21900 

25200 

T static (K) 

370 

430 

510 

Mach no. 

3.5 

4.0 

4.4 


5.3 Modeling Approach 

Simplifying assumptions must be made in order to model a gas-particle flow through the facility 
mixer, the facility nozzle, along two external compression ramps, through the inlet throat area and the 
isolator section using a quasi- ID code. Since the facility mixer and nozzle are assumed to be 
axisymmetric, it is straightforward to cast their geometry into a quasi- ID format. However, the 
precompression plate, final external inlet ramp and isolator geometries are inherently two-dimensional, as 
are the oblique shocks preceding each of those components and the expansion fan leading into the 
isolator. This difficulty was addressed by modeling each of these last three components separately from 
the mixer/nozzle section. The jump conditions across each of the oblique shocks were calculated from the 
exit conditions of the upstream component and the specified turning angle using classical methods, such 
as those found in Saad (1985). The calculated post-shock conditions were then used as the inflow 
boundary conditions for the next component. For the isolator, the post-shock conditions were used as the 
input to a classical expansion fan calculation to determine the inflow boundary conditions of that 
component. Each of these last components was then treated as a constant area section having a length 
equal to the corresponding body-side bounding surface, as shown in Figure 31. The isolator was 
considered as a constant area section, even though a slight area increase with length is shown in 
Figure 26. However, that area relief is designed to accommodate boundary layer growth so that the 
isolator can behave as a truly constant area section and not prematurely choke near the end of the duct. 
Since the 1-D analysis does not include boundary layer, a constant area duct was modeled. The cross- 
sectional area for these constant area sections can be any arbitrary value, though a value of 0.1 nr was 
generally used throughout. It should be noted that the length of each component is representative of the 
gas-particle flow near the bounding wall. Gas-particle flows further out from the wall will have a shorter 
residence time between shocks. 

One additional modeling step was required for the isolator. Isolator flow is generally characterized by 
a complex, three-dimensional series of reflected oblique shocks, possibly followed by a weakened normal 
shock in the ramjet operating mode. It is not realistic to try and model this flowfield using a 1-D code. 
Instead, a uniform negative momentum source term was added to the momentum Equation (2), at each 
grid point along the isolator length to reflect the total pressure loss generated by the internal shock 
structure. This approach yielded nearly linear changes in the primary flow variables between the isolator 
inlet and exit, ending at the desired isolator exit values given in Table 3. This behavior for the various 
flow parameters can be observed in the various plots in the later Results section. 

With these modeling assumptions, a 1 percent reacting particle mass fraction case had a typical total 
run time of approximately 80 hr on a dual-processor (2 GHz) PC workstation. This relatively long run 
time was driven primarily by the low convective velocity in the low-speed mixer section. A similar case 
with just a single reacting particle typically took approximately 40 CPU hours on the same machine. In 
both cases, a converged no-particle solution was used as the initial condition. 
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Figure 32. — Nozzle area distribution where X = 0 indicates the beginning of 
the convergent section of the nozzle. 

5.4 Grid Descriptions 

Simple uniform 1-D grids were used for each component. For the mixer/nozzle simulation, an overall 
grid totaling 4000 points was split two sections. An 800 point uniform grid was used for the 300 cm 
constant-area mixer section and a 3200 point uniform grid was used for the 300 cm nozzle section. The 
exact area distribution used for all of the simulations for the 3 m long nozzle, as shown schematically in 
Figure 26, is given in Figure 32. Each of the remaining three components downstream of the nozzle was 
modeled using uniform 300 point grids. 

6.0 Results 

6.1 Inert Particle Calculations 

A series of Q1D calculations were performed with inert alumina particles utilizing the test facility and 
scramjet geometries described in the previous section at the Mach 5 and 6 total temperature flow 
conditions defined therein. The Mach 7 total temperature condition was not run because this condition is 
near the melting point of alumina (2323 K). As mentioned in the introduction, particle diameters of 10 pm 
to 1 mm (1000 pm) were considered for this study. In addition to examining the two extreme particle 
diameters, an intermediate diameter of 1 00 pm was also included in the study to help define the effect of 
particle diameter on particle and gas flow conditions entering the model scramjet combustor. The 
majority of results shown below are for single particles traversing the facility mixer/nozzle, 
precompression plate, final inlet compression ramp and isolator components. Additionally, a single case 
with a 1 percent particle mass fraction was run to determine if the effect of increased particle mass 
fraction on the scramjet combustor inflow conditions warranted additional calculations at elevated inert 
particle mass fractions. 

The particles used in this portion of the study were assumed to be low porosity, high purity tabular 
alumina. All the required material properties, including temperature dependent properties, were available 
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from Munro (1997), with the exception of particle emissivity. However, Pluchino and Masturzo (1981) 
determined that the emissivity of small alumina particles approaches zero as the material purity 
approaches 1 00 percent. As the tabular alumina is a very high purity form of alumina, it was therefore 
assumed that the particle emissivity was in fact zero for these calculations, eliminating all radiative heat 
transfer from the particles to the water-cooled facility walls (500 K wall temperature). The published 
density for tabular alumina is 3984 Kg/nT. Temperature dependant values of heat capacity and thermal 
conductivity are shown in Figures 33 and 34, respectively. 




Temperature (K) 

Figure 34. — Tabular alumina thermal conductivity. 
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It should be noted that these property values had to be extrapolated above 1800 K. However, as each 
of these properties approach an asymptote by that temperature, no significant error is incurred. For 
particles with an internal temperature gradient, mass-averaged values of these temperature dependent 
properties were calculated and applied at each time step. 

Prior to running these inert particle cases, it was desirable to determine whether or not it would be 
necessary to include gas-phase chemistry in the calculations. Gas-only calculations with and without gas- 
phase chemistry were run at the Mach 5 and 6 total temperature conditions. The Mach 5 and 6 
calculations showed no significant difference in the bulk fluid properties between the results with and 
without gas-phase chemistry. All minor species in the finite-rate chemistry calculations remained below 1 
part-per-million throughout the computational domain for both conditions. Based on these results, it was 
determined that the inert particle cases could reasonably be run without using finite rate chemistry. By 
running all of the inert particle cases in this manner, significant computational run-time improvements 
were achieved. 

A total of six cases were run using single 1 0-, 1 00- and 1 000-pm alumina particles at the Mach 5 and 
6 total temperature flow conditions. Initial test runs showed that the Biot number for the 10 pm particles 
never exceeded 0. 1 during the simulation, indicating that those particles should be modeled as isothermal. 
The larger particles were modeled with internal grids of 5 1 grid points for the 100 pm particles and 101 
grid points for the 1000 pm particles. The number of internal grid points for each particle size was 
selected by running a series of test cases with increasingly refined grids until the solution became grid 
independent. Figures 35 to 38 show the gas and particle surface temperature and velocity results for these 
six cases, each plot at a single facility operating point total temperature. Only a single gas-phase solution 
is given in each of these figures as the gas-phase solutions at any given Mach number are not effected by 
the insertion of a single particle. The discontinuities in the gas flow results indicate the positions of the 
oblique shock waves indicated in Figure 3 1 at each flow turning point. The particle initial temperature 
was set 200 K higher than the initial gas temperature in the mixer section to simulate the overheating of 
the test gas typically required in storage heater facility operations, as described earlier. The particle initial 
velocity is set to the incoming gas flow velocity. 

Examining these figures, several observations can be made. First, that while the 10 pm particles 
reasonably follow the gas flow temperature and velocity profiles, the 100 and 1000 pm particles 
significantly lag the gas flow due to their higher mass to surface area ratio that decreases the particle drag 
relative to the particle inertia. Second, that while the 1 0 pm particles respond to the step-change in gas 
flow properties at the oblique shock positions, albeit with significant delay, the 1 00 and 1 000 pm particles 
show almost no response to the shocks for the same reason as the first observation above. And lastly, that 
increasing the facility total temperature from the Mach 5 condition to the Mach 6 condition does little to 
alter the response of the particles, other than to shift the temperature and velocity scales. This shift in 
scales is more evident in Figures 39 to 44, where the data is replotted keeping the particle size constant in 
each figure, while varying the simulated Mach number total temperature. 
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Precompression Plate 



Axial Distance (m) 

Figure 35. — Comparison of gas and single particle surface temperature profiles for 
various particle sizes for the Mach 5 total temperature flow condition. 



Figure 36. — Comparison of gas and single particle velocity profiles for 
various particle sizes for the Mach 5 total temperature flow condition. 
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Figure 37. — Comparison of gas and single particle surface temperature profiles 
for various particle sizes for the Mach 6 total temperature flow condition. 



Figure 38. — Comparison of gas and single particle velocity profiles for 
various particle sizes for the Mach 6 total temperature flow condition. 
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Figure 39. — Comparison of gas and 10 pm particle surface temperature 
profiles for Mach 5 and 6 total temperature flow conditions. 



Figure 40. — Comparison of gas and 100 pm particle surface temperature 
profiles for Mach 5 and 6 total temperature flow conditions. 
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Figure 41 . — Comparison of gas and 1000 pm particle surface temperature 
profiles for Mach 5 and 6 total temperature flow conditions. 



Figure 42. — Comparison of gas and 10 pm particle velocity profiles for 
Mach 5 and 6 total temperature flow conditions. 


N ASA/TM— 20 10-216765 


51 


2000 


o 

o 

•c 

> 


1750- 
1500 - 
1250 - 
1000 - 
750 - 
500 - 
250 - 



V -Mach 5 

V -Mach 5 

V -Mach 6 

V -Mach 6 

p 


T 


T 


Axial Distance (m) 

Figure 43. — Comparison of gas and 100 pm particle velocity profiles for 
Mach 5 and 6 total temperature flow conditions. 



Mach 5 and 6 total temperature flow conditions. 
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Figures 45 and 46 show the internal particle temperature profiles for the 100 pm particles and 
1 000 pm particles, respectively, at the exit plane of the mixer-nozzle section (MN) and the exit plan of 
the isolator (IS). The temperature values are normalized by temperature at which they are initially inserted 
in the flow stream, T imert , while the particle radius coordinate, r, is normalized by the particle outer radius, 
r p . It is of interest that the 100 pm particle internal temperature profiles shift somewhat with the change in 
simulated Mach number in Figure 45, indicative of significant particle temperature change, but that the 
particle internal temperature profiles remain almost constant with the change in simulated Mach number 
for the 1000 pm particles shown in Figure 46. It is also evident, as expected, that the larger particles, with 
their larger thermal mass, sustain a larger thermal gradient at the particle surface. 

Lastly, the effect of a relatively high particle mass fraction is examined in Figures 47 to 52. In 
Figure 47, the gas and particle temperature profiles for a case using a single 10 pm particle are compared 
to the gas and particle temperature profiles for a case wherein the flow contained a 1 percent mass 
fraction of 10 pm particles inserted into the flow at the same conditions as the single particles. The 
profiles are nearly identical for both the gas and particles. Figures 48 and 49 show a closer view of the 
nozzle exit and the ramp and isolator sections where a small difference in gas and particle temperatures is 
evident, however. After the same manner, the overall gas and particle velocity profiles are shown in 
Figure 50, with more detailed views given in Figures 51 and 52. Looking at these results, one must 
conclude that while some difference does exist in gas and particle properties due to the presence of a 
1 percent mass fraction of particles, the differences are insignificant. Based on this result, no further 
1 percent inert particle mass fraction simulations were performed. 
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Figure 45. — Radial temperature profiles of 100 pm particles for the Mach 5 
and 6 conditions at the exit planes of the mixer/nozzle (MN) and the 
isolator (IS). 
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Figure 46. — Radial temperature profiles of 1000 pm particles for the Mach 5 
and 6 conditions at the exit planes of the mixer/nozzle (MN) and the 
isolator (IS). 
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Figure 47. — Comparison of gas and particle temperature profiles for the 
Mach 6 total temperature flow condition fora single 10 pm particle and 
for a 1 percent 10 pm particle mass fraction. 
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Figure 48. — Comparison of gas temperature profiles for the Mach 6 total 
temperature flow condition for a single 1 0 pm particle and for a 1 
percent particle mass fraction of 1 0 pm particles over the last 1 .0 m of 
the nozzle, the precompression plate, the final external compression 
ramp, and the isolator. 
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total temperature flow condition for a single 10 pm particle and for a 
1 percent particle mass fraction of 10 pm particles over the last 1.0 m of 
the nozzle, the precompression plate, the final external compression 
ramp, and the isolator. 
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Figure 50. — Comparison of gas and particle velocity profiles for the Mach 
6 condition fora single 10 pm particle and fora 1 percent 10 pm particle 
mass fraction. 



Figure 51 . — Comparison of gas velocity profiles for the Mach 6 total 
temperature flow condition for a single 1 0 pm particle and for a 
1 percent particle mass fraction of 10 pm particles over the last 1.0 m 
of the nozzle, the precompression plate, the final external compression 
ramp, and the isolator. 
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Figure 52. — Comparison of particle velocity profiles for the Mach 6 total 
temperature flow condition for a single 1 0 pm particle and for a 
1 percent particle mass fraction of 10 pm particles over the last 1.0 m 
of the nozzle, the precompression plate, the final external compression 
ramp, and the isolator. 

6.2 Single Graphite Particle Calculations 

Using the kinetic and physical property information for graphite previously described, the Mach 5 and 
6 total temperature single inert particle simulations were repeated using single graphite particles of the 
same particle diameters. For graphite the radiative emissivity is non-zero, so radiation heat transfer from 
the hot particles to the facility wall is included in each graphite particle simulation. A value for the 
emissivity of graphite of 0.9 was used (Pierson, 1993). Also, since the allowable temperature for graphite 
is considerably higher than that for alumina, Mach 7 total temperature simulations were performed for all 
three particle sizes. As with the inert particle cases, gas-phase finite rate chemistry was not used for the 
Mach 5 and 6 calculations. For the Mach 7 total temperature simulations, examination of gas-phase only 
calculations with and without finite-rate chemistry showed that significant mass fractions of minor 
species, particularly NO, were generated, as shown in Figure 53. 

The formation of NO in turn had a measurable impact upon the 0 2 and N 2 concentrations, as shown in 
Figure 54, which in turn had an effect on the static temperature, shown in Figure 55. Based on these results, 
the Mach 7 total temperature simulations were run with gas-phase finite-rate chemistry. In order to decrease 
computational time, those reactions given in Appendix B that include hydrogen-containing species were 
eliminated for these calculations, as there is no hydrogen in the system up through the isolator section. 

One additional adjustment was made for the graphite particle calculations in comparison to the alumina 
particle calculations. For the alumina particles, an internal grid was used for both the 100 and 1000 pm 
diameter particles. Flowever, the much higher thermal conductivity of graphite compared to alumina 
resulted in a Biot number for the 1 00 pm graphite particle below 0. 1 during the entire particle trajectory. 
During the code validation portion of this study, it was shown that such low Biot numbers, particularly at 
high heat flux conditions, yielded lower accuracy heat transfer results. Therefore, the 100 pm graphite 
particle cases were run without an internal particle grid (i.e. with the whole particle at a single temperature). 
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Figure 53. — Mass fractions of minor species from Mach 7 total 
temperature reacting flow simulation 
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Figure 54. — Comparison of mass fractions of major species from Mach 7 total 
temperature simulations with and without finite-rate chemistry. 
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Figure 55. — Comparison of static temperature profiles from Mach 7 total 
temperature simulations with and without finite-rate chemistry. 
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Figures 56 and 57 show the particle temperature results for the Mach 5 and Mach 6 total temperature 
conditions, respectively, along with the corresponding inert alumina particle results for comparison. 

Figure 58 shows the temperature results for the Mach 7 condition that was not run for the inert particles. 
The 1 0 pm graphite particle is immediately quenched at the Mach 5 condition by the colder gas around it 
due to its high surface area to mass ratio leading to a high heat loss to the surroundings relative to the 
combustion heat addition at this condition, and stays quenched throughout the particle trajectory. The 
different thermal response compared to the 10 pm alumina particle in the nozzle is due to the differences 
in density, porosity, heat capacity, and thermal conductivity between the two materials. For the Mach 6 
and 7 simulations, the 10 pm particle is consumed within a very short distance as the combustion heat 
addition is now much greater than the heat loss to the environment. The code does not contain a specific 
particle burnout model, and so the final phases of particle combustion may be somewhat inaccurate, but 
considering how rapidly the particle bums out, it is unlikely that a more accurate modeling approach 
would yield a substantially different result. For all cases, burnout is assumed to occur when the particle 
density reaches 1 percent of the initial particle density, at which point the particle is removed from the 
simulation. 

The 100 pm graphite particles show similar behavior for the Mach 5, 6 and 7 total temperature 
conditions. In the constant-area low speed mixer, the first three meters, combustion is diffusion limited, as 
the particle and gas velocities are nearly identical. However, as the particle enters the convergent portion 
of the nozzle, a particle velocity lag is introduced, which greatly increases the transport of oxygen to the 
particle surface, and the temperature rapidly increases. The temperature rise then levels off, and for the 
Mach 6 case even drops off somewhat, based on the competing effects of lower free stream pressure in 
the nozzle leading to lower oxygen concentration, decreased particle mass due to combustion, increased 
particle internal surface area, heat transfer and combustion kinetics. For the Mach 7 total temperature 
condition, the particle burned out just after entering the convergent portion of the nozzle. 
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Figure 56. — Mach 5 graphite versus alumina surface temperature comparison. 



Figure 57. — Mach 6 graphite versus alumina surface temperature comparison. 


N ASA/TM— 20 10-216765 


60 



The 1000 pm particles show very little temperature change at any of the three simulated Mach 
numbers due to their size. The relatively low external surface area to particle mass ratio of the 1000 pm 
particles decreases the impact of radiation and convective heat transfer on particle temperature. 
Simultaneously, intra-particle diffusion of oxygen and combustion products is greatly restricted, limiting 
combustion to the outer regions of the particle. 

Figure 59 shows the internal particle temperature profiles for the 1000 pm particle at the exit plane of 
the mixer-nozzle section (MN) and the exit plan of the isolator (IS). The T/T inseit values above 1.0 are 
indicative of ongoing particle combustion. 

Figures 60 to 62 show the velocity results for all three simulated Mach numbers and all three particle 
sizes, with the Mach 5 and Mach 6 results compared to their corresponding alumina particle cases. The 
Mach 6 and Mach 7 10 pm particle results are not readily discernable in these plots as the particles bum 
out so quickly that their trajectories do not differentiate themselves from the other cases. Since graphite is 
significantly less dense than alumina, it is not surprising to see that the graphite particles more closely 
follow the gas flow than the corresponding alumina particles do. This effect is accentuated due to mass 
loss from oxidation of the graphite particles, particularly for the 100 pm particles. In Figure 63, graphite 
particle porosity histories are given for each particle size and simulated Mach number. This plot shows 
that at each simulated Mach number, the 100 pm particles approach complete burnout, allowing the 
100 pm graphite particles to respond much more quickly to changes in flow velocity than the 
corresponding alumina particles. The 1000 pm graphite particles lose a relatively small amount of mass 
traversing the domain in each case, and so behave similarly to the alumina particles, but reach a higher 
velocity due to the density difference. 
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Figure 59. — Radial temperature profiles of the 1 000 pm particle for the 
Mach 5, 6 and 7 conditions at the exit planes of the mixer/nozzle (MN) 
and the isolator (IS). 
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Figure 60. — Mach 5 graphite versus alumina velocity comparison. 
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Figure 61. — Mach 6 graphite versus alumina velocity comparison. 



It should be noted that for the 1 00 and 1 000 pm graphite particles, the increase in porosity, 
corresponding to mass loss, does not increase as much as might be expected with simulated Mach 
number, despite the higher gas and particle temperatures. This is due to two causes. First, the particle 
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combustion is primarily diffusion limited throughout the particle trajectory. The higher temperature cases 
do not have significantly higher oxygen diffusion rates to and within the particle, so the mass loss cannot 
proceed at a much higher rate despite the increase in kinetic rate with temperature. Second, an increase in 
radiative heat transfer from the particle surface with increasing temperature tends to offset the internal 
heat generation from increased combustion rate. 

To gain additional insight into the various factors effecting particle velocity and temperature profiles 
as a function of particle size, it is of interest to examine some additional parameters. Figure 64 give 
particle-gas relative Mach number (defined as V gas -V partic i e divided by the local sound speed) data for all 
three particle sizes at the Mach 5 total temperature condition. The Mach 5 data was selected here because 
at the higher Mach numbers the 1 0 pm particle quickly burned out and therefore these values could not be 
compared throughout the particle trajectories. The relative Mach number traces shown in Figure 64 trend 
as expected, with the larger, heavier particles showing a larger velocity deficit. The 1 000 pm particles 
have a supersonic relative Mach number through much of the nozzle, while the smaller particles never 
achieve a supersonic relative Mach number. While Mach number is generally defined as a scalar, for 
these Mach number calculations the relative Mach number is allowed to go negative when the particle 
velocity is greater than the gas flow velocity. Therefore, negative Mach numbers are shown in Figure 64, 
occurring in each case for the 10 pm particle after a shock jump-condition has been imposed on the gas 
flow and for the last two segments of the flowpath for the 1 00 pm particle. 

Figures 65 to 67 show the different heat flows into, out of, and accumulating in the particle outer 
shell, as defined in Figure 3, for the three Mach 5 total temperature conditions. For the 10 and 100 pm 
particles, the “outer shell” actually encompasses the entire particle, since no internal grid was used. This 
means that the conduction term into or out of the outer shell is by definition zero, and is not shown. For 
the 1 000 pm particle, a 1 0 1 point internal grid was used, so the outer shell has a thickness equal to 1 
percent of the particle radius (5 pm). It should be noted that in these figures, some of the heat flow terms 
have been multiplied by 1 0 to make them visible on the scale of the plot. 


1 - 

o.; 

0.8 

0.7 

0.6 

0.5 

0.4 

0.3 



Mach 5 

Mach 6 
Mach 7 
Mach 5 
Mach 6 
Mach 7 
Mach 5 
Mach 6 
Mach 7 


1 0 Micron 

10 Micron 
1 0 Micron 
100 Micron 
100 Micron 
100 Micron 
1000 Micron 
1 000 Micron 
1000 Micron 




0.2 


Axial Distance (m) 

Figure 63. — Graphite particle porosity. 
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Figure 64. — Mach 5 gas-particle relative mach number. 


Looking at Figure 65, heat flows for the 10 pm particle, we see high initial heat flow rates as the 
particle cools from the particle insertion temperature (200 K higher than the gas temperature) down to the 
gas temperature. The particle quickly reaches steady state within the mixer section, up to the 3.0 m point, 
and the heat flows stay constant until the particle enters the nozzle section. As the gas temperature drops 
in the nozzle, combustion is quenched and the particle radiation drops to zero. At the same time, 
convection spikes, removing heat from the particle as shown by the negative heat accumulation term. 
Finally, there is some particle heating across each of the oblique shocks and throughout the isolator 
section, located at approximately 6.0 m, 7.0 m, and 8.0 to 8.4 m, respectively. 

For the 100 pm particle, Figure 66, we see that the combustion term is always larger than the 
convection heat loss term (convection out of the particle is defined as being positive, as is radiation), 
resulting in continuous particle temperature rise, as the radiation heat loss remains small. In the area of 
the isolator, the sum of the convection and radiation heat flows do approach the value of the combustion 
term, so the temperature remains fairly constant, as shown in Figure 56. 
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Figure 65. — Heat flows for 10 pm graphite particle at Mach 5 total temperature condition. 
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Figure 66. — Heat flows for 100 pm graphite particle at Mach 5 total 
temperature condition. 
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Figure 67. — Heat flows for 1000 gm graphite particle outer shell at Mach 5 
total temperature condition. 


Lastly, the heat flows (including internal conduction) for the 1 000 pm graphite particle at the Mach 5 
condition are given in Figure 67. As this particle’s temperature changes very little throughout the 
simulation, we see that the radiative heat loss is essentially constant. As with the other particle sizes, the 
convection term increases rapidly entering the nozzle, but the actual change in particle temperature is 
small because of the large thermal mass represented by this relatively large particle. The ramp shocks and 
isolator section are discemable in Figure 67, but have only a very minor effect on the particle temperature 
profile. 

6.3 One Percent Mass Fraction Graphite Particle Simulations 

The same simulations that were run for the inert particles and the single graphite particles were then 
run using a 1 percent mass fraction of graphite particles. All of the 1 percent mass fraction simulations 
were run using finite rate chemistry in the gas phase. As in the Mach 7 single graphite particle 
simulations, the hydrogen-containing species were eliminated from the kinetic mechanism to speed the 
calculations. Each 1 percent mass fraction simulation was run until a time invariant solution was obtained 
for both the particle and gas streams. 

Figures 68 to 76 show the temperature and velocity histories of 10, 100, and 1000 pm graphite 
particles and the corresponding gas stream solutions calculated using a 1 percent mass fraction of 
particles at each of the three simulated Mach numbers. These figures also contain the single graphite 
particle results for comparison. 

For the 10 pm particles at the Mach 5 total temperature flow condition, Figure 68, there is virtually no 
difference between the single particle case and the 1 percent mass fraction case. Since combustion of the 
particle is almost immediately quenched by the cooler surrounding gas, there is no impact on the gas 
stream that could cause a change in the overall solution. The insertion into the gas stream of a 1 percent 
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ass fraction of particles at a temperature 200 K higher than the gas stream is insufficient to cause a 
significant change in the gas stream variables, as was observed with the inert alumina particles. At the 
Mach 6 and Mach 7 total temperature flow conditions, Figures 69 and 70, the 10 pm graphite particles 
undergo very rapid combustion. The heat release from the particles significantly raises the gas stream 
temperature in the mixer. This increase in mixer temperature results in an increased isolator exit 
temperature, as well as increased velocity in sections of the flow path. The increase in velocity is due to 
an increase in the sound speed due to the higher static temperature in regions of approximately constant 
Mach number. The somewhat lower relative temperature increase in the Mach 7 case, Figure 70, is due 
primarily to the increase in gas specific heat with increasing temperature. 



single 10 pm graphite particle and a 1 percent mass fraction of 10 pm graphite 
particles at the Mach 5 total temperature flow condition. 
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Figure 69. — Comparison of particle temperature and velocity profiles 

between a single 10 pm graphite particle and a 1 percent mass fraction of 
10 pm graphite particles at the Mach 6 total temperature flow condition. 
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Figure 70. — Comparison of particle temperature and velocity profiles 

between a single 10 pm graphite particle and a 1 percent mass fraction of 
10 pm graphite particles at the Mach 7 total temperature flow condition. 
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Comparing the results for the Mach 5 and Mach 6 total temperature flow condition for the 1 00 pm 
particles, Figures 71 and 72, similar results are observed. While the single particles do not bum out within 
the computational domain, the higher gas temperature for the 1 percent mass fraction simulations from 
graphite combustion leads to more rapid particle consumption and eventual burnout. As would be 
expected, the particles bum out more rapidly at the higher Mach 6 total temperature condition. 

Some unique features are visible in the Mach 7 results for the 1 00 pm graphite particles, as shown in 
Figure 73. Flere the single graphite particle and the 1 percent graphite particle stream bum out at 
approximately the same point. This occurs despite the increasing gas temperature found in the 1 percent 
mass fraction results. This occurs because at this very high temperature condition the combustion process 
is completely diffusion limited. Increasing the gas temperature does not translate into a higher particle 
combustion rate. 

Lastly, in Figures 74 to 76, little difference is observed between the single particle solutions and the 
1 percent mass fraction solutions for the 1 000 pm particles at each of the three simulated Mach numbers. 
Since relatively little combustion is occurring in these cases, the gas streams are unaffected by the 
presence of the 1 percent mass fraction of particles. 

Figures 77 to 79 show the particle porosity histories for the single graphite particles and the 1 percent 
mass fraction graphite particle stream as a function of particle size and simulated Mach number. 



Figure 71 . — Comparison of particle temperature and velocity profiles between 
a single 100 pm graphite particle and a 1 percent mass fraction of 100 pm 
graphite particles at the Mach 5 total temperature flow condition. 
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Figure 72. — Comparison of particle temperature and velocity profiles between 
a single 100 pm graphite particle and a 1 percent mass fraction of 100 pm 
graphite particles at the Mach 6 total temperature flow condition. 
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Figure 73. — Comparison of particle temperature and velocity profiles between 
a single 100 pm graphite particle and a 1 percent mass fraction of 100 pm 
graphite particles at the Mach 7 total temperature flow condition. 
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Figure 74. — Comparison of particle temperature and velocity profiles between 
a single 1000 pm graphite particle and a 1 percent mass fraction of 
1000 pm graphite particles at the Mach 5 total temperature flow condition. 



Figure 75. — Comparison of particle temperature and velocity profiles between 
a single 1000 pm graphite particle and a 1 percent mass fraction of 
1000 pm graphite particles at the Mach 6 total temperature flow condition. 
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Figure 76. — Comparison of particle temperature and velocity profiles between 
a single 1000 pm graphite particle and a 1 percent mass fraction of 
1000 pm graphite particles at the Mach 7 total temperature flow condition. 



Figure 77. — Comparison of particle porosity profiles between single graphite 
particles and a 1 percent mass fraction of graphite particles at the Mach 5 
total temperature flow condition. 
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Figure 78. — Comparison of particle porosity profiles between single graphite 
particles and a 1 percent mass fraction of graphite particles at the Mach 6 
total temperature flow condition. 
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Figure 79. — Comparison of particle porosity profiles between single graphite 
particles and a 1 percent mass fraction of graphite particles at the Mach 7 
total temperature flow condition. 
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Figure 80 shows the NO mass fraction profile for all the 1 percent graphite particle mass fraction 
cases where NO is produced above a mass fraction of lxlO -6 at any point in the computational domain. 
All of the Mach 5 total temperature flow condition cases and the Mach 6 total temperature flow condition 
case with no particles have NO levels below this threshold. The Mach 6 cases with particles all have very 
low NO mass fractions which only make it above the lx 10 6 mass fraction threshold because of particle 
combustion. However, all of the Mach 7 simulations have significant mass fractions of NO. The NO is 
generated by the high temperature in the mixing section, particularly when enhanced by particle 
combustion. The NO level decreases as the temperature drops with expansion, but the NO reactions are 
eventually “frozen” in the expansion section of the nozzle due to the rapid temperature drop along with 
relatively slow nitrogen chemistry. Looking at Figures 81 and 82, the N0 2 and O mass fraction profiles, it 
is evident that a small mass fraction of O is generated by the high temperatures in the mixing section in 
addition to the NO, and that this atomic oxygen combines with the NO to produce N0 2 during the 
expansion just upstream and downstream of the nozzle throat. These N0 2 mass fractions are then also 
“frozen” as they progress downstream. Prominent in Figure 81 is a spike in the N0 2 mass fraction early in 
the mixer for the 1 00 pm graphite particles at the Mach 7 total temperature flow condition. This spike 
represents a region wherein the kinetics favors the formation of N0 2 from the NO and O formed at the 
initial high temperature. As the temperature increases due to graphite combustion, the kinetics transitions 
to favor the formation of NO and O instead of N0 2 , and the N0 2 previously formed dissociates back to 
NO and O. In general, it can be concluded from these figures that the more rapid the particle combustion 
process, the higher the levels of NO and O that are generated, and therefore an increased mass fraction of 
N0 2 is generated in the expansion. 



Axial Distance (m) 

Figure 80. — Comparison of NO mass fraction profiles between no graphite particles 
and a 1 percent mass fraction of graphite particles. 
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Figure 81 . — Comparison of NO 2 mass fraction profiles between no graphite 
particles and a 1 percent mass fraction of graphite particles. 
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Figure 82. — Comparison of O mass fraction profiles between no graphite 
particles and a 1 percent mass fraction of graphite particles. 
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Figures 83 and 84 show the spatial distribution of C0 2 and CO, respectively. As expected, the 
combinations of particle size and simulated Mach number that yield little or no particle combustion show 
low levels of C0 2 and CO. This would include all the 1000 pm particle cases and the Mach 5, 10 pm 
particle case. Low temperature graphite combustion favors CO production at the particle surface and 
inhibits CO oxidation in the gas phase, as seen with the Mach 5 total temperature, 100 pm graphite 
particle case. As the temperature increases due to more rapid graphite combustion, such as with the Mach 
6 total temperature, 10 and 100 pm particle cases, CO is initially produced by the particle, but this 
transitions to C0 2 production with rising particle temperature. With increasing gas phase temperature, the 
CO is oxidized to C0 2 . This gas phase CO oxidation takes much longer for the slower burning Mach 6, 
100 pm particle case, and is not completed before the nozzle expansion lowers the flow temperature. At 
that point, CO production is again favored over C0 2 at the particle surface, and the CO level rises through 
the nozzle. This rise in the CO mass fraction for the Mach 6, 100 pm case in the nozzle is also contributed 
to by a small amount of C0 2 dissociation, as indicated by the drop in C0 2 concentration in Figure 83. 
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Figure 83. — Comparison of C0 2 mass fraction profiles between 1 percent 
graphite particle mass fraction cases. 
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Figure 84. — Comparison of CO mass fraction profiles between 1 percent 
graphite particle mass fraction cases. 

6.4 Ignition Calculations 

As a final step in examining the effect of graphite particulate vitiation on hypersonic propulsion 
testing, a series of ignition delay calculations were performed based on the isolator section exit flow 
conditions for each 1 percent graphite particle mass fraction case. Hydrogen was used as the fuel at an 
equivalence ratio of 0.8 based upon the amount of molecular oxygen present at the exit of the isolator. 

The computational domain generally consisted of a constant area duct 0.5 meters long, open at each end. 
450 equally spaced grid points were used in each model. In some cases the length of the domain was 
extended to 1.0 m if ignition did not occur within the 0.5 m length. A supersonic inlet boundary condition 
was used with the pressure, temperature, velocity, and species mass fractions taken from the isolator exit 
conditions, as given in Tables 6 and 7. If particles were present at the isolator exit, they were included in 
the ignition delay calculation at their isolator exit condition, Table 8, and were allowed to continue to 
react based on the local flow conditions. The exit boundary conditions remained supersonic for each case. 
The initial condition within the duct was 100 percent N 2 at the temperature, pressure, and velocity of the 
isolator exit. Pure nitrogen gas was used as the initial condition for the calculation to help distinguish the 
location of the interface between the gas coming from the in-flow boundary versus and the gas initially in 
the duct for this plug-flow simulation. Ignition was defined as the point in time when a 200 K temperature 
rise above the inlet boundary condition was observed within the computational domain. Figure 85 shows 
a typical result at the point of ignition. 


TABLE 6.— GAS TEMPERATURES (K) AT THE EXIT OF THE ISOLATOR SECTION 


Mach no. 

No particles 

10 pm 

100 pm 

1000 pm 

5 

798 

800 

908 

802 

6 

1008 

1229 

1201 

1022 

7 

1264 

1452 

1441 

1278 
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TABLE 7.— SPECIES MASS FRACTIONS AT THE EXIT OF THE ISOLATOR SECTION 


Mach no./ 
particle size 

o 2 

n 2 

NO 

no 2 

n 

o 

CO 

O 


0.233000 

0.767000 













0.233000 

0.767000 












7 / None 

0.229035 

0.763448 

0.007340 

0.000177 




5/10 

0.232141 

0.766549 



0.001175 

0.000135 


6/10 

0.207284 

0.757233 

0.000347 

0.000018 

0.035111 

0.000007 


7/10 

0.196415 

0.743837 

0.019493 

0.000323 

0.035401 

0.000141 


5/100 

0.218665 

0.756865 



0.002132 

0.022322 

0.000016 

6/100 

0.209855 

0.756875 

0.000084 

0.000008 

0.026443 

0.006705 


7/100 

0.196023 

0.744401 

0.017977 

0.000304 

0.035747 

0.000668 


5 / 1000 

0.232698 

0.766849 



0.000271 

0.000182 


6/1000 

0.232200 

0.766511 

0.000002 


0.000359 

0.000923 

0.000005 

7/1000 

0.227262 

0.762568 

0.007643 

0.000183 

0.001472 

0.000872 



TABLE 8.— PARTICLE SURFACE TEMPERATURE, VELOCITY, AND POROSITY 


FOR SURVIVING GRAPHITE PARTICLES AT THE EXIT OF THE ISOLATOR SECTION 


Mach no./ 
particle size 

Surface temperature 
(K) 

Velocity 

(m/s) 

Porosity 

5/10 

625 

1115 

0.264 

5 / 1000 

1387 

785 

0.247 

6 / 1000 

1911 

941 

0.273 

7/1000 

2421 

1156 

0.294 



Axial Distance (rn) 

Figure 85. — Example of ignition point for a IVvitiated air/particle mixture from the 
Mach 6 total temperature simulation with a 1 percent mass fraction of 1000 pm 
graphite particles. 
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The calculated ignition delay times for each condition with a 1 percent initial mass fraction of 
graphite particles, along with the ignition delay times for simulations done without any particles, are 
given in Table 9. In those cases where little graphite combustion occurred, only a slight reduction in 
ignition is observed, or even a slight increase in ignition delay (due to the heat-sink effect of the particles 
during the ignition delay calculation) for the Mach 5, 10 pm particle case. However, in the calculations 
where significant graphite combustion has occurred, the ignition delay time is greatly reduced from the 
cases with no particles. To examine this effect further, a second set of ignition delay calculations were 
performed to determine if the decrease in ignition delay time was due to vitiate species being present 
(including particles) or graphite combustion heat release. These calculations were performed at the 
elevated temperatures given in Table 6, but using no gas phase or particulate vitiates, only regular air. 
These results are given in Table 10. Comparing the results from Tables 9 and 10, it is evident that the 
majority of the graphite particle vitiation effect for these conditions is simply heat release. The presence 
of the vitiate species in the flow generally serves to slightly retard ignition. The relatively small impact of 
the combined vitiate species is consistent with other studies, such as Han, et al. (2004), where the 
significant effects of known ignition enhancing species such as NO are generally seen at lower 
temperatures and higher pressures than those studied here, and where the expected effects of a particular 
vitiate is frequently cancelled out by the effects of other vitiate species when they are studied in 
combination. The purpose of these ignition delay calculations is only to look at the combined effect of all 
the vitiates. Further studies of individual vitiate species effects would be needed to determine more 
detailed information, but such studies are not warranted here. 


TABLE 9.— IGNITION DELAY TIME IN MICROSECONDS WITH VITIATES 


Mach no. 

No particles 

10 pm 

100 pm 

1000 pm 

5 

1105 

1112 

367 

1055 

6 

183 

82.7 

91.2 

169 

7 

56.1 

43.1 

41.9 

52.1 

TABLE 10.— 

IGNITION DELAY TIME IN MICROSECONDS WITHOUT VITIATES 

Mach no. 

No particles 

10 pm 

100 pm 

1000 pm 

5 

1105 

1094 

353 

1055 

6 

183 

75.8 

84.4 

169 

7 

56.0 

39.1 

37.8 

51.7 


7.0 Conclusions 

In order to improve the understanding of particulate vitiation in hypersonic propulsion test facilities, a 
quasi-one dimensional numerical tool has been developed capable of efficiently modeling reacting 
particle-gas flows over a wide range of gas and particle conditions. Features of this code include finite 
rate kinetics for gas phase chemistry, a global porous-particle combustion model for the solid phase, 
mass, momentum and energy interactions between phases, subsonic and supersonic particle drag and heat 
transfer models, variable time step, and appropriate sets of subsonic and supersonic boundary conditions. 
The basic capabilities of this tool have been validated against available literature data or other validated 
codes. 

This code has been demonstrated through a parametric study of particle vitiation effects in a model 
hypersonic propulsion test facility and scramjet. Parameters studied were simulated flight Mach number 
(Mach 5, 6, and 7), particle size (10, 100, and 1000 pm diameters), particle mass fraction (single particles 
and 0.01) and particle material (alumina and graphite). For the inert alumina particles, it was found that 
the presence of particles up to the 1 percent mass fraction had very little effect on the gas phase. As a 1 
percent mass fraction would in reality be an extreme case, it is safe to conclude that hot wall ignition on 
the particle surface is likely to be the only potential inert particle vitiation effect, an effect that was 
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beyond the scope of this study. It was observed that only the 1 0 pm alumina particles tracked the gas flow 
reasonably well, and that when the particle size was increased to 1 000 pm, the particle responded very 
little at all to changes in the gas flow, including shock waves, even though supersonic relative Mach 
numbers are observed. This is to be expected as the large particles have significant inertia against velocity 
changes, and have a relatively large volume -to-surface-area ratio that decreases heat transfer effects with 
the surroundings. 

With the graphite particles, it was observed that the smallest particles in the study will either quickly 
quench, or quickly be consumed, depending on the gas and particle temperatures. The combustion of 
these 1 0 pm particles is solely determined by the particle temperature -dependent kinetic rates and not by 
oxygen diffusion. As the particle size increases, the particle is less likely to quench, because it has a larger 
thermal inertia, and so will remain at a high temperature long enough for the particle to ignite. As the 
particle and surrounding gas temperatures increase, the particle consumption rate increases somewhat, but 
not as fast as might be expected. This is due to diffusion rate limitations, both to the particle surface and 
internal to the particle. For the largest particles, combustion is significantly diffusion limited, so particle 
and gas temperature have little effect on the combustion rate. As the particle mass fraction is stepped up 
to 1 percent, the main change is the addition of significant heat release. As the choice of 1 percent mass 
fraction was meant to be a limiting case, the amount of heat release is extreme and not expected to be 
observed in real facilities. However, it does serve to increase the particle combustion rate, and to increase 
the quantity of NOx generated in the flow. In those cases where low reaction rates of graphite were 
observed for single particles, the increase to 1 percent mass fraction had very little impact, similar to the 
inert particles. Ignition delay time calculations for the 1 percent mass fraction of graphite particles cases 
showed significant decreases in ignition delay in the cases where higher graphite combustion rates, and 
therefore higher gas temperatures, were observed, but further calculations showed that this was due 
primarily to the increased combustor inlet temperature, not the vitiate species present in the flow. What 
chemical effects were present appeared to slightly increase the ignition delay. A review of several 
summary papers on gas-phase vitiate effects did little to clarify this result, as the findings are typically 
highly case-specific and not easily compared. 

While the results discussed above are by no means conclusive regarding the potential effects of 
particle vitiation, they do provide significant insight into the physical effects of gas flow conditions on the 
particles, and of the particles on the gas flow. Understanding of the relative effects of the particle drag, 
the various forms of heat transfer, and combustion throughout the test facility and test engine is an 
important first step. The tool developed as part of this effort can be useful in determining the likely effects 
of specific particle contamination situations at specific test conditions, when the required modeling inputs 
are made available. This numerical tool is also generally useful for modeling any gas phase or gas-solid 
two-phase fluid flow problem that can be modeled as quasi-one dimensional. 

8.0 Recommendations for Future Work 

During the course of demonstrating the numerical tool developed for this study, a number of potential 
improvements were noted. These can be broadly categorized into three areas, as follows. 

8.1 Modeling Improvements 

1 . Particle burnout model. — A particle burnout model should be added to the code if medium size 
particles are to be modeled. This is not important for very small particles which bum out quickly or for 
very large particles which suffer low mass loss. 

2. Particle diameter shrinkage model. — This is another model that would be most critical for medium 
size particles. For smaller particles that bum out or quench quickly, the impact of incorrectly modeling 
the particle shrinkage is small, as there would not be expected to be a significant difference in the 
numerical solution at any significant distance downstream of the insertion point. For large particles, little 


N ASA/TM— 20 10-216765 


81 



diameter change would be expected unless extremely long residence times were being considered. A 
particle diameter shrinkage model would need to be based on empirical data for the material, particle 
geometry, and test conditions under study. 

3. Multiple simultaneous particle sizes. — The code currently only models a single particle size for a 
given run. The capability for handling multiple particle sizes simultaneously would enhance the ability of 
the code to model real situations. 

4. Hot wall ignition model. — The ability to examine gaseous fuel-oxidizer ignition at the surface of an 
individual hot, reacting particle could be important. This type of modeling is most likely beyond the 
capability of a quasi-one dimensional code. A multi-dimensional code that addresses multi-species 
diffusion to and from the particle surface and particle surface chemistry would be needed to complement 
the current capability. 

8.2 Input Data Improvements 

1. Graphite Combustion Data. — For accurate results with reacting particles, it is necessary to have 
graphite combustion data, including combustion rate and particle shrinkage, across the full range of 
conditions to be experienced by the particles in the model. Such data could then be applied to the global 
reaction model in the code. 

2. Good particle characterization. — For modeling a real system, it is necessary to have characterized 
the particle contamination in the facility. Particle size distribution and particle mass fraction must be 
known. Also the particles themselves need to be fully characterized in terms of overall porosity and 
internal pore size. 

3. Validation Data. — ft was not possible to validate the quasi-one dimensional code as a whole 
because of a lack of relevant validation data. Each capability within the code was validated separately, to 
the greatest extent possible. Gathering applicable validation data would greatly enhance confidence in the 
accuracy of the numerical results. 

8.3 Additional Modeling 

1. Lower Particle Mass Fractions for Reacting Particles. — The upper limit of 1 percent particle mass 
fraction used for this study was unrealistically high, although it does provide an upper bound. Cases could 
be run at various particle mass fractions below this level to understand the effect of particle vitiation at 
more realistic levels for those facility test conditions at which particle vitiation effects were observed. 

2. Additional Particle Sizes for Reacting Particles. — In this demonstration study, different modes of 
behavior were observed at the different particle sizes. Running more studies in the intermediate size range 
would be of significant interest, since it is in this size range that the particle behavior is the least 
predictable. 
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Appendix A. — Symbols List 

Gas phase mixture specific heat 
Specific heat of gaseous species i 
Particle internal pore diameter 
Gas phase specific internal energy 
Effectiveness factor coefficient 
Gas phase Gibbs energy of formation 

Gibbs energy of formation of gaseous species i 

Gas phase specific enthalpy 

Specific enthalpy of gaseous species i 

Particle convection coefficient 

Gas phase enthalpy of formation 

Enthalpy of formation of gaseous species i 

Backward reaction rate 

Forward reaction rate 

Gas thermal conductivity 

1 st kinetic rate for carbon combustion 

Observed carbon combustion reaction order 

Arrhenius reaction temperature exponent 

Number of chemical species 

Pressure 

Heat accumulation within the particle outer shell 
Gas phase conductive heat flux 

Conductive heat flow into the particle from the particle outer shell 

Convective heat flow out of the particle outer shell 

Combustion heat addition in particle outer shell 

Radiative heat flow out of the particle outer shell 

Wall heat flux 

Particle radial coordinate 

Particle combustion rate based on particle outer shell energy balance 
Mass-averaged particle combustion rate based on carbon kinetics 
Local internal particle combustion rate based on carbon kinetics 
Particle combustion rate based on oxygen diffusion 
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r p Particle radius 

r* Inner product used to calculate artificial viscosity 

t Time 

u g Gas phase velocity 

u p Particle velocity 

x Axial distance coordinate 

x p Particle axial position 

A Cross sectional area 

Af Arrhenius reaction pre-exponential factor 

C Artificial viscosity interim variable 

Cj) Particle drag coefficient 

C s Particle material heat capacity 

D e Particle internal diffusion coefficient 

Dj Diffusion coefficient of gaseous species i 

Dj j Binary diffusion coefficient between gaseous species i and j 

D k Knudsen diffusion coefficient 

D ( -) 2 m Molecular diffusion coefficient of oxygen into gas mixture 

E a Arrhenius reaction activation energy 

F Array of convective terms from the gas phase governing equations 

F p Drag force on a particle 

H\ Scramjet inlet throat height 

J Array of source terms from the gas phase governing equations 

K d Mass transfer coefficient 

K eq Equilibrium constant 

K P 2 nd kinetic rate for carbon combustion 

K s Particle thermal conductivity 

K ± Artificial viscosity interim variable 

Li Scramjet inlet and isolator ramp lengths 

M Mixture average molecular weight 

M c Carbon atomic mass 

Mi Species molecular weight 

MR Ratio of moles of CO to moles of CO 2 produced in carbon combustion 

M re i Mach number based on particle/gas relative velocity 

MSR Molecular speed ratio 

M t Thiele modulus 
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N Number of chemical reactions 

Nurf Nusselt number based on particle diameter 

Pq 2 f s Free stream partial pressure of oxygen 

P() l S Partial pressure of oxygen at the particle surface 

Pr Prandtl number 

R Mixture gas constant 

Re,. e / Particle Reynolds number based on particle/gas relative velocity 
R u Universal gas constant 

S e Energy equation energy source term 

S ext Particle external surface area 

Sh Sherwood number 

Si Reactant or product species 

Sint Particle internal surface area 

S m Continuity equation mass source term 

S m j Species continuity equation mass source term 

S p Momentum equation momentum source term 

Stot Particle total surface area 

Tfn m Particle film temperature 

T g Gas temperature 

T g *j j Reduced gas temperature for collision between species i and j 

T insert The temperature at which a particle is inserted into the flow stream 

T p Particle surface temperature or bulk temperature for uniform temperature particles 

T p i Particle temperature at internal grid point i 

T rec Convective recovery temperature 

Twall Facility mixer/nozzle wall temperature 

T t: j Effective temperature of gaseous species i 

T s i j Average effective temperature between gaseous species i and j 

T *2,2 Reduced temperature for mixture viscosity 

U Array of dependent variables from the gas phase governing equations 

Xj Mole fraction of gaseous species i 

Yj Mass fraction of gaseous species i 

a Ratio of particle mass flow to gas mass flow 

a Scramjet inlet/isolator ramp lengths 
P,- Scramjet inlet/isolator shock angles 
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Artificial viscosity 

Particle porosity 

Particle emissivity 

Gas phase volume fraction 

Averaging variable for mixture viscosity 

Ratio of specific heats 

Particle internal diffusion effectiveness factor 
Mass-averaged particle internal diffusion effectiveness factor 

Artificial viscosity velocity eigenvalue 
Gas phase viscosity 
Viscosity of gaseous species i 
Courant number 

Diffusion velocity of gaseous species i 
Stoichiometric coefficient of reactant species i 
Stoichiometric coefficient of product species i 

Stoichiometric coefficient of oxygen consumed per mole of carbon consumed 
Artificial viscosity flux limiter 
Gas phase density 

Particle material density 

Effective collision diameter of gaseous species i 

Average effective collision diameter between species i and j 

Stefan-Boltzmann constant 

Oxygen diffusion coefficient 

Numerical time step 

Particle internal grid size 

Collision integral between species i and j 

Collision integral for mixture viscosity 

Intermediate variable for carbon combustion rate calculation 

Ratio of moles of CO2 to total moles produced in carbon combustion 
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Appendix B. — Gas Phase Kinetic Mechanism 

The detailed gas phase kinetic mechanism used for this study is given in Table B. 1. While the 
majority of the reactions follow the simple Arrhenius form, there are a number that have additional 
features. Ten of the reactions required the participation of a third body, for which the species third body 
efficiencies are included in Table B.2. Reactions 18, 19, 23, and 25 include a second duplicate reaction 
with different Arrhenius constants. The sum of the two reactions provides the correct temperature 
dependent behavior. Lastly, reactions 12, 14, and 24 are pressure dependent reactions. The first set of 
Arrhenius constants represents the high pressure limit, while the second set represents the low pressure 
limit. These two reaction rates are then combined via the Troe method, using the constants provided, for 
reactions 12 and 14, and via the Lindemann method for reaction 24. These methods are described in detail 
by Kee, et al (1993). 


TABLE B.I.— GAS-PHASE KINETIC MECHANISM" 


Reaction 
ref. no. 

Reaction 

A f 

n 

E a 

Reaction 

source 3 

1 

h+o 2 ^O+OH 

2.644xl0 16 

-0.6707 

17041 

05Davis 

2 

o+h 2 <=>h+oh 

4.589xl0 4 

2.70 

6260 

05Davis 

3 

oh+h 2 ^h+h 2 o 

1.734xl0 8 

1.51 

3430 

05Davis 

4 

20H <=> 0 + H 2 0 

3.973xl0 4 

2.40 

-2110 

05Davis 

5 

H + H + M <=> H2 + M 

1.78xl0 18 

-1.0 

0.0 

05 Davis 

6 

2H + H 2 <=> 2H 2 

9.0xl0 16 

-0.600 

0.0 

05 Davis 

7 

2H + H 2 0<=>H 2 +H 2 0 

5.624xl0 19 

-1.25 

0.0 

05Davis 

8 

2H + C0 2 <=>H2 + C0 2 

5.50xl0 20 

-2.0 

0.0 

05Davis 

9 

h+oh+m<=>h 2 o+m 

4.40xl0 22 

-2.0 

0.0 

05 Davis 

10 

O + H + M <=> OH + M 

9.428xl0 18 

-1.0 

0.0 

05Davis 

11 

O + O + M <=> 0 2 +M 

1.20xl0 17 

-1.0 

0.0 

05Davis 

12 

H + 0 2 (+M) <» H0 2 (+M) 

5.116xl0 12 

0.44 

0.0 

05 Davis 


Low 

6.328xl0 19 

-1.4 

0.0 



TROE 

0.5 

lxlO -30 

lxlO 30 


13 

h 2 + o 2 <=> ho 2 +h 

5.916x10 s 

2.433 

53502 

05Davis 

14 5 

20H(+m)<=> H 2 0 2 (+m) 

1 . 1 1 x 10 14 

-0.37 

0.0 

05 Davis 


Low 

2.01xl0 17 

-0.584 

-2293 



TROE 

0.0 

345 

10 


15 

ho 2 + H <=> 0 + h 2 o 

3.97xl0 12 

0.0 

671 

05Davis 

16 

H0 2 +H <=> 20H 

7.485xl0 13 

0.0 

295 

05Davis 

17 

HO 2 + O <^> OH + O 2 

4xl0 13 

0.0 

0.0 

05 Davis 

18a 

ho 2 + oh<=>o 2 +h 2 o 

2.375xl0 13 

0.0 

-500 

05Davis 

18b 

ho 2 +oh»o 2 +h 2 o 

lxlO 16 

0.0 

17330 

05Davis 

19a 

2H02 02 +H 202 

1.3x10“ 

0.0 

-1630 

05Davis 

19b 

2H02 < w > 02 +H 202 

3.658xl0 14 

0.0 

12000 

05Davis 

20 

h 2 o 2 + h <=> ho 2 + h 2 

6.05xl0 6 

2.0 

5200 

05 Davis 

21 

h 2 o 2 +h»oh+h 2 o 

2.41xl0 13 

0.0 

3970 

05 Davis 

22 

h 2 o 2 + 0»0L1 + H0 2 

9.63xl0 6 

2.0 

3970 

05Davis 

23a 

h 2 o 2 +oh-»ho 2 +h 2 o 

2xl0 12 

0.0 

427 

05Davis 

23b 

h 2 o 2 +oh-»ho 2 +h 2 o 

2.67xl0 41 

-7.0 

37600 

05Davis 

24 

co+o(+m)<=>co 2 (+m) 

1.362xl0 10 

0.0 

2384 

05Davis 


Low 

1.173xl0 24 

-2.79 

4191 


25a 

C0+0HOC02+H 

8x10“ 

0.14 

7352 

05Davis 

25b 

CO + OH 0 CO 2 + H 

8.784xl0 10 

0.03 

-16 

05Davis 
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TABLE B.I.— GAS-PHASE KINETIC MECHANISM 0 


Reaction 
ref. no. 

Reaction 

A f 

n 

E a 

Reaction 

source 3 

26 

co+o 2 «co 2 +o 

1.119xl0 12 

0.0 

47700 

05Davis 

27 

CO + HO, «C0 2 +OH 

3.01x10° 

0.0 

23000 

05Davis 

28 

hco+h»co+h 2 

1.2xl0 14 

0.0 

0.0 

05Davis 

29 

HCO + O <» CO + OH 

3x10° 

0 

0 

05Davis 

30 

HCO + OoC0 2 +II 

3x10° 

0 

0 

05Davis 

31 

hco+oh<=>co+h 2 o 

3.02x10° 

0 

0 

05Davis 

32 

HCO + M<=>CO + H + M 

1.87x10° 

-1 

17000 

05Davis 

33 

HCO + H 2 0 <=> CO + H + H20 

2.244x10° 

-1 

17000 

05Davis 

34 

hco+o 2 <=>co+ho 2 

1.204xl0 10 

.807 

-727 

05Davis 

35 

n + no«n 2 +0 

2.7x10° 

0 

355 

GRI3.0 

36 

N + 0 2 »N0 + 0 

9xl0 9 

1 

6500 

GRI3.0 

37 

N +OH <=> NO + H 

3.36x10° 

0 

385 

GRI3.0 

38 

H0 2 +N0«N0 2 +0H 

2.11x10° 

0 

-480 

GRI3.0 

39 

NO + O + M <=> N0 2 +M 

1.06xl0 20 

-1.4 

0 

GRI3.0 

40 

N0 2 +0 <=> N0 + 0 2 

3.9x10° 

0 

-240 

GRI3.0 

41 

no 2 +h<=>no+oh 

1.32xl0 14 

0 

360 

GRI3.0 

42 

H + NO + M <=> HNO + M 

4.48xl0 19 

-1.3 

740 

GRI3.0 

43 

HNO + O <=> NO + OH 

2.5x10° 

0 

0 

GRI3.0 

44 

HNO+Hotl, +NO 

9x10“ 

.7 

660 

GRI3.0 

45 

HNO + OHoNO + H 2 0 

1.3xl0 7 

1.9 

-950 

GRI3.0 

46 

hno+o 2 »ho 2 +no 

1x10° 

0 

13000 

GRI3.0 

47 

N + C0 2 »NO + CO 

3x10° 

0 

11300 

GRI3.0 


a 05Davis refers to the mechanism developed by Davis, et al (2005), while GRI3.0 refers to the GRI 3.0 (2005) mechanism. 
b For this reaction, the optional fourth Troe constant is used with a value of 345. 

c The units for the pre-exponential factor Aj arc Kg/m 2 -s-atm, while the units for the activation energy E„ are cal/g-mole. The temperature 
exponent n is unit-less. 


TABLE B.2.— THIRD BODY EFFICIENCIES 


Species 

React 

5 

React 

9 

React 

10 

React 

11 

React 

12 

React 

14 

React 

24 

React 

32 

React 

39 

React 

42 

h 2 

0 

2 

2 

2.4 

.75 

2 

2 

2 

1 

1 

o 2 

1 

1 

1 

1 

.85 

1 

1 

1 

1 

1 

H,0 

0 

6.3 

12 

15.4 

11.89 

6 

12 

0 

1 

1 

OH 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

H 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

O 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

ho 2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

h 2 o 2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

N 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

NO 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

no 2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

HNO 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

co 2 

0 

3.6 

3.6 

3.6 

2.18 

3.6 

3.6 

3.6 

1 

1 

CO 

1 

1.75 

1.75 

1.75 

1.09 

1.75 

1.75 

1.75 

1 

1 

HCO 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

Ar 

.63 

.38 

.7 

.83 

.4 

.7 

.7 

1 

1 

1 

n 2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 
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